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1.  Introduction 


1.1  Description  of  Problem 

Dynamic  failure  in  bound  particulate  materials  is  a  combination  of  physical  processes 
including  grain  and  matrix  deformation,  intra-granular  cracking,  matrix  cracking,  and 
inter- granular- matrix/binder  cracking/debonding,  and  is  influenced  by  global  initial 
boundary  value  problem  (IBVP)  conditions.  Discovering  how  these  processes  occur  by 
experimental  measurements  is  difficult  because  of  their  dynamic  nature  and  the  influence  of 
global  boundary  conditions  (BCs).  Typically,  post-mortem  microscopy  observations  are 
made  of  fractured/fragmented/comminuted  material  (Kipp  et  ah,  1993),  or  real-time 
in-situ  infrared-optical  surface  observations  are  conducted  of  the  dynamic  failure  process 
(Guduru  et  ah,  2001).  These  observation  techniques,  however,  miss  the  origins  of  dynamic 
failure  internally  in  the  material.  Under  quasi-static  loading  conditions,  non- destructive 
high  spatial  resolution  (a  few  microns)  synchrotron  micro-computed  tomography  can  be 
conducted  (Fredrich  et  ah,  2006)*  to  track  three-dimensionally  the  internal  grain-scale 
fracture  process  leading  to  macro-cracks  (though  these  cracks  can  propagate  unstably). 
Dynamic  loading,  however,  can  generate  a  significantly  different  microstructural  response, 
usually  fragmented  and  comminuted  material  (Kipp  et  ah,  1993).  Global  BCs,  such  as 
lateral  confinement  on  cylindrical  compression  specimens,  also  can  influence  the  resulting 
failure  mode,  generating  in  a  glass  ceramic  composite  axial  splitting  and  fragmentation 
when  there  is  no  confinement  and  shear  fractures  with  confinement  (Chen  and 
Ravichandran,  1997).  Thus,  we  resort  to  physics-based  modeling  to  help  uncover  these 
origins  dynamically. 

Examples  of  bound  particulate  materials  include,  but  are  not  limited  to,  the  following: 
polycrystalline  ceramics  (crystalline  grains  with  amorphous  grain  boundary  phases,  figure 
1(a)),  metal  matrix  composites  (metallic  grains  with  bulk  amorphous  metallic  binder,  figure 
1(b)),  particulate  energetic  materials  (explosive  crystalline  grains  with  polymeric  binder, 
figure  1(c)),  asphalt  pavement  (stone/rubber  aggregate  with  hardened  binder,  figure  1(d)), 
mortar  (sand  grains  with  cement  binder),  conventional  quasi-brittle  concrete  (stone 
aggregate  with  cement  binder),  and  sandstones  (sand  grains  with  clayey  binder).  Bound 
particulate  materials  contain  grains  (quasi-brittle  or  ductile)  bound  by  binder  material 
often  called  the  “matrix.”  The  heterogeneous  particulate  nature  of  these  materials  governs 
their  mechanical  behavior  at  the  grain-to-macroscales,  especially  in  IBVPs  for  which 
localized  deformation  nucleates.  Thus,  grain-scale  material  model  resolution  is  needed  in 
regions  of  localized  deformation  nucleation  (e.g.,  at  a  macro-crack  tip,  or  at  the  high  shear 

*Such  experimental  techniques  are  not  yet  mature,  but  can  provide  meaningful  insight  into  the  origins  of 
‘static’  fracture,  and  thus  could  play  an  important  role  in  the  discovery  of  the  origins  of  dynamic  failure. 
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(c)  (d) 


Figure  1.  (a)  Microstructure  of  alumina,  composed  of  grains  bound  by  a  glassy  phase. 

(b)  Silicon  carbide  (SiC)  reinforced  2080  aluminum  metal  matrix  composite 
(Chawla,  2004).  The  four  black  squares  are  indents  to  identify  the  region,  (c) 
Cracking  in  explosive  HMX  grains  and  at  the  grain-matrix  interfaces  (Baer, 
2007).  (d)  Cracking  in  asphalt  pavement. 


strain  rate  interface  region  between  a  projectile  and  target  material').  To  predict  dynamic 
failure  for  realistic  IBVPs,  a  modeling  approach  needs  to  account  simultaneously  for  the 
underlying  grain-scale  physics  and  macroscale  continuum  IBVP  conditions. 

Traditional  single-scale  continuum  constitutive  models  have  provided  the  basis  for 
understanding  the  dynamic  failure  of  these  materials  for  IBVPs  on  the  macroscale 
(Rajendran  and  Grove,  1996;  Dienes  et  ah,  2006;  Johnson  and  Holmquist,  1999),  but 
cannot  predict  dynamic  failure  because  they  do  not  account  explicitly  for  the  material’s 
particulate  nature.  Direct  Numerical  Simulation  (DNS)  directly  represents  the  grain-scale 
mechanical  behavior  under  static  (Caballero  et  ah,  2006)  and  dynamic  loading  conditions 
(Kraft  et  ah,  2008;  Kraft  and  Molinari,  2008b).  Currently,  DNS  is  the  best  approach  to 
understanding  fundamentally  dynamic  material  failure,  but  IT  is  deficient  in  the  following 
ways:  (1)  it  is  limited  by  current  computing  power  (even  massively  parallel  computing)  to 
a  small  representative  volume  element  (RVE)  of  the  material  and  (2)  it  usually  must 
assume  unrealistic  BCs  on  the  RVE  (e.g.,  periodic,  or  prescribed  uniform  traction  or 

^Both  projectile  and  target  material  could  be  modeled  with  such  grain-scale  material  model  resolution 
at  their  interface  region  where  significant  fracture  and  comminution  occurs.  We  start  by  assuming  the 
projectile  is  a  deformable  solid  continuum  body  without  grain-scale  resolution,  and  then  extend  to  include 
such  resolution  in  the  future. 
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displacement).  Thus,  multiscale  modeling  techniques  are  needed  to  predict  dynamic  failure 
in  bound  particulate  materials. 

Current  multiscale  approaches  attempt  to  do  this,  but  fall  short  by  one  or  more  of  the 
following  limitations:  (1)  not  providing  proper  BCs  on  the  microstructural  DNS  region 
(called  the  “unit  cell”  by  Feyel  and  Chaboche  (2000),  extended  to  account  for 
discontinuities  in  Bclytschko  et  al.  (2008));  (2)  homogenizing  at  the  macroscale  the 
underlying  microstructural  response  in  the  unit  cell  and  thus  not  maintaining  a 
computational  “open  window”  to  model  microstructurally  dynamic  failure*;  and  (3)  not 
making  these  methods  adaptive,  i.e.,  moving  a  computational  “open  window”  with 
grain-scale  model  resolution  over  regions  experiencing  dynamic  failure. 

Feyel  and  Chaboche  (2000)  and  Belytschko  et  al.  (2008)  recognized  the  complexities  and 
limitations  of  unit  cell  methods  as  they  are  currently  formulated,  implemented,  and  applied. 
Feyel  (2003)  stated  that,  in  addition  to  the  periodicity  assumption  for  the  microstructure 
(impossible  to  model  fracture),  “...  real  structures  have  edges,  either  external  or  internal 
ones  (in  case  of  a  multimaterial  structure) .  In  the  present  FE2  framework,  nothing  has 
been  done  to  treat  such  effects.  As  a  consequence,  one  cannot  expect  a  good  solution  near 
edges.  This  is  clearly  a  weak  point  of  the  approach  ...”  In  fact,  for  a  non-periodic 
heterogeneous  microstructure  found  in  bound  particulate  materials,  we  should  not  expect 
predictive  results  for  modeling  nucleation  of  fracture  anywhere  in  the  unit  cell. 

Bclytschko  et  al.  (2008)  introduced  discontinuities  into  Feyel  and  Chaboche’s  (2000)  unit 
cell  (calling  it  a  “perforated  unit  cell”)  and  relaxed  the  periodicity  assumption  to  model 
fracture  nucleation,  while  up-scaling  the  effects  of  unit  cell  discontinuities  to  the  macroscale 
to  obtain  global  cracks  embedded  in  the  finite  element  (FE)  solution  (using  the  extended 
FE  method).  BCs  on  the  unit  cell  are  an  issue,  as  well  as  the  interaction  of  adjacent  unit 
cells.  As  noted  in  Belytschko  et  al.  (2008),  if  regular  displacement  BCs  (i.e.,  no  jumps)  are 
applied  to  unit  cells  that  are  fracturing,  then  the  fracture  is  constrained  non-physically. 
Bclytschko  et  al.  (2008)  proposed  to  address  this  issue  by  solving  iteratively  for 
displacement  BCs  by  applying  a  traction  instead.  What  traction  to  apply  is  still  an 
unknown  and  can  be  provided  by  the  coarse-scale  FE  solution.  Belytschko  et  al.  (2008) 
stated  that  “...  the  application  of  boundary  conditions  on  the  unit  cell  and  information 
transfer  to/from  the  unit  cell  pose  several  difficulties  ...  When  the  unit  cell  localizes, 
prescribed  linear  displacements  as  given  in  the  analysis  are  not  compatible  with  the 
discontinuities  ...  The  effects  of  boundaries  and  adjacent  discontinuities  are  not  reflected  in 
the  method.” 

Whis  is  a  problem  especially  for  modeling  fragmentation  and  comminution  microstructurally. 
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1.2  Proposed  Approach 


A  finite  strain  micromorphic  plasticity  model  framework  (Regueiro,  2010b)  is  applied  to 
formulate  a  simple  pressure-sensitive  plasticity  model  to  account  for  the  underlying 
microstructural  mechanical  response  in  bound  particulate  materials  (pressure-sensitive 
heterogeneous  materials).  Linear  isotropic  elasticity  and  nonassociative  Drucker-Prager 
plasticity  with  cohesion  hardening/softening  are  assumed  for  the  constitutive  equations 
(Regueiro,  2009).  Micromorphic  continuum  mechanics  is  used  in  the  sense  of  Eringen 
(1999).  This  was  found  to  be  one  of  the  more  general  higher  order  continuum  mechanics 
frameworks  for  accounting  for  underlying  microstructural  mechanical  response.  Until  this 
work,  however,  the  finite  strain  formulation  based  on  multiplicative  decomposition  of  the 
deformation  gradient  F  and  microdeformation  tensor  y  has  not  been  presented  in  the 
literature  with  sufficient  account  of  the  reduced  dissipation  inequality  and  conjugate  plastic 
power  terms  to  dictate  the  plastic  evolution  equation  forms.  We  provide  such  details  in  this 
report. 

To  illustrate  the  application  of  the  micromorphic  plasticity  model  to  the  problem  of 
interest,  we  refer  figure  2,  which  illustrates  a  concurrent  multiscale  modeling  framework  for 
bound  particulate  materials  (target)  impacted  by  a  deformable  solid  (projectile).  The 
higher  order  continuum  micromorphic  plasticity  model  is  used  in  the  overlap  region 
between  a  continuum  FE  and  DNS  representation  of  the  particulate  material.  The 
additional  degrees  of  freedom  provided  by  the  micromorphic  model  (microshear, 
micro-dilation/compaction,  and  microrotation)  allow  the  overlap  region  to  be  placed  closer 
to  the  region  of  interest,  such  as  at  a  projectile-target  interface.  Further,  from  this 
interface  region,  standard  continuum  mechanics  and  constitutive  models  can  be  used. 

1.3  Focus  of  Report 

Regarding  the  approach  described  in  section  1.2,  this  report  primarily  on  the  nonlinear 
micromorphic  continuum  mechanics  and  finite  strain  elastoplasticity  constitutive  model 
tasks.  How  this  generalized  continuum  model  couples  via  an  overlapping  region  to  the  DNS 
region  (figure  2)  is  described  in  sections  2.4  and  2.5. 

The  discrete  element  (DE)  and/or  FE  representation  of  the  particulate  microstructure  is 
intentionally  not  shown  so  as  not  to  clutter  the  drawing  of  the  microstructure.  The  grains 
(binder  matrix  not  shown)  of  the  microstructure  are  “meshed”  using  DEs  and/or  FEs  with 
cohesive  surface  elements  (CSEs).  The  open  circles  denote  continuum  FE  nodes  that  have 
prescribed  degrees  of  freedom  (DOFs)  D  based  on  the  underlying  grain-scale  response, 
while  the  solid  circles  denote  continuum  FE  nodes  that  have  free  DOFs  D  governed  by  the 
micromorphic  continuum  model.  We  intentionally  leave  an  “open  window”  (i.e.,  DNS)  on 
the  particulate  microstructural  mesh  in  order  to  model  dynamic  failure.  If  the  continuum 
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Figure  2.  Two-dimensional  illustration  of  the  concurrent  computational  multiscale  mod¬ 
eling  approach  in  the  contact  interface  region  between  a  bound  particulate  ma¬ 
terial  (e.g.,  ceramic  target)  and  deformable  solid  body  (e.g.,  refractory  metal 
projectile). 

mesh  overlays  the  whole  particulate  microstructural  region,  as  in  Klein  and  Zimmerman 
(2006)  for  atomistic-continuum  coupling,  then  the  continuum  FEs  would  eventually  become 
too  deformed  by  following  the  microstructural  motion  during  fragmentation.  The 
blue-dashed  box  at  the  bottom-center  of  the  illustration  is  a  micromorphic  continuum  FE 
region  that  can  be  converted  to  a  DNS  region  for  adaptive  high-fidelity  material  modeling 
as  the  projectile  penetrates  the  target. 

An  outline  of  the  report  is  as  follows:  section  2.1  summarizes  the  statement  of  work  (SOW) 
and  the  tasks,  2.2  presents  the  formulation  of  the  nonlinear  (finite  deformation) 
micromorphic  continuum  mechanics  balance  equations,  sedetion  2.3  presents  the  finite 
strain  elastoplasticity  modeling  framework  based  on  a  multiplicative  decomposition  of  the 
deformation  gradient  and  microdeformation  tensor,  sections  2.4  and  2.5  describe  how  the 
micromorphic  continuum  mechanics  fits  into  a  multiscale  modeling  approach,  section  3 
summarizes  the  results,  and  section  4  provides  the  conclusions  and  future  work. 

1.4  Notation 

Cartesian  coordinates  are  assumed  for  easier  presentation  of  concepts  and  also  in  order  to 
define  a  Lagrangian  elastic  strain  measure  Ee  in  the  intermediate  configuration  B , 
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assuming  a  multiplicative  decomposition  of  the  deformation  gradient  F  and 
microdeformation  tensor  x  into  elastic  and  plastic  parts  (section  2.3.1).  See  Regueiro 
(2010)  for  more  details  regarding  finite  strain  micromorphic  elastoplasticity  in  general 
curvilinear  coordinates,  and  also  Eringen  (1962)  for  nonlinear  continuum  mechanics  in 
general  curvilinear  coordinates  and  Clayton  et  al.  (2004,  2005)  for  nonlinear  crystal 
elastoplasticity  in  general  curvilinear  coordinates. 

Index  notation  is  used  so  as  to  be  as  clear  as  possible  with  regard  to  details  of  the 
formulation.  Cartesian  coordinates  are  assumed,  so  all  indices  are  subscripts,  and  spatial 
partial  derivative  is  the  same  as  covariant  derivative  (Eringen,  1962).  Some  symbolic/direct 
notation  is  also  given,  such  that  (ab)ifc  =  u, (a  0  h)ijki  =  ai:jbki,  (a  0  c)ijk  =  aimcjmk. 
Boldface  denotes  a  tensor  or  vector,  where  its  index  notation  is  given.  Generally,  variables 
in  uppercase  letters  and  no  overbar  live  in  the  reference  configuration  B$  (such  as  the 
reference  differential  volume  dV),  variables  in  lowercase  live  in  the  current  configuration  B 
(such  as  the  current  differential  volume  dv) ,  and  variables  in  uppercase  with  overbar  live  in 
the  intermediate  configuration  B  (such  as  the  intermediate  differential  volume  dV).  The 
same  applies  to  their  indices,  such  that  a  differential  line  segment  in  the  current 
configuration  dxi  is  related  to  a  differential  line  segment  in  the  reference  configuration  dXj 
through  the  deformation  gradient:  dxi  =  FudXj  (Einstein’s  summation  convention 
assumed  [see  Eringen,  1962;  Holzapfel,  2000]).  In  addition,  the  multiplicative 
decomposition  of  the  deformation  gradient  is  written  as  Fik  =  F^Fjr  (F  =  Ff  Fp),  where 
superscripts  e  and  p  denote  clastic  and  plastic  parts,  respectively.  Subscripts  (•),;  (•)  /  and 
(•),/  imply  spatial  partial  derivatives  in  the  current,  intermediate,  and  reference 
configurations,  respectively.  A  superscript  prime  symbol  (•)'  denotes  a  variable  associated 
with  the  microelement  for  micromorphic  continuum  mechanics.  Superposed  dot 
(□)  =  D(\3)/Dt  denotes  material  time  derivative.  The  symbol  =f  implies  a  definition. 


2.  Technical  Discussion 


2.1  Statement  of  Work  (SOW)  and  Specific  Tasks 

Bound  particulate  materials  are  commonly  found  in  industrial  products,  construction 
materials,  and  nature  (e.g.,  geological  materials).  They  include  polycrystalline  ceramics 
(e.g.,  crystalline  grains  with  amorphous  grain  boundary  phases),  energetic  materials  (high 
explosives  and  solid  rocket  propellant),  hot  asphalt,  asphalt  pavement  (after  asphalt  has 
cured),  mortar,  conventional  quasi- brittle  concrete,  ductile  fiber  composite  concretes,  and 
sandstones,  for  instance.  Bound  particulate  materials  contain  particles^  (quasi-brittle  or 
ductile)  bound  by  binder  material  often  called  the  “matrix.” 

§We  use  “particle”  and  “grain”  interchangeably. 
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The  heterogeneous  nature  of  bound  particulate  materials  governs  its  mechanical  behavior 
at  the  particle-  to  continuum- scales.  The  particle-scale  is  denoted  as  the  scale  at  which 
particle-matrix  mechanical  behavior  is  dominant,  thus  necessitating  that  particles  and 
matrix  material  be  resolved  explicitly  (i.e.,  meshed  directly  in  a  numerical  model), 
accounting  for  their  interfaces  and  differences  in  material  properties.  Currently,  there  is  no 
approach  enabling  prediction  of  initiation  and  propagation  of  dynamic  fracture  in  bound 
particulate  materials — for  example,  polycrystalline  ceramics,  particulate  energetic 
materials,  mortar,  and  sandstone — accounting  for  their  underlying  particulate 
microstructure  across  multiple  length-scales  concurrently.  Traditional  continuum  methods 
have  provided  the  basis  for  understanding  the  dynamic  fracture  of  these  materials,  but 
cannot  predict  the  initiation  of  dynamic  fracture  without  accounting  for  the  material’s 
particulate  nature.  DNS  of  deformation,  intra-particle  cracking,  and 
inter-particle- matrix /binder  debonding  at  the  particle-scale  is  limited  by  current 
computing  power  (even  massively  parallel  computing)  to  a  small  RVE  of  the  material,  and 
usually  must  assume  overly- restrictive  BCs  on  the  RVE  (e.g.,  fixed  normal  displacement). 

Multiscale  modeling  techniques  are  clearly  needed  to  accurately  capture  the  response  of 
bound  particulate  materials  in  a  way  accounting  simultaneously  for  effects  of  the 
microstructure  at  the  particle-scale  and  BCs  applied  to  the  engineering  structure  of 
interest,  at  the  continuum-scale.  The  services  of  a  scientist  or  engineer  are  required  to 
develop  the  mathematical  theory  and  numerical  methodology  for  multiscale  modeling  of 
bound  particulate  materials  of  interest  to  the  U.S.  Army  Research  Laboratory  (ARL). 

The  overall  objective  of  the  proposed  research  is  to  develop  a  concurrent  multiscale 
computational  modeling  approach  that  couples  regions  of  continuum  deformation  to 
regions  of  particle-matrix  deformation,  cracking,  and  debonding,  while  bridging  the 
particle-  to  continuum-scale  mechanics  to  allow  numerical  adaptivity  in  modeling  initiation 
of  dynamic  fracture  and  degradation  in  bound  particulate  materials. 

For  computational  efficiency,  the  solicited  research  use  DNSs  only  in  the  spatial  regions  of 
interest,  such  as  the  initiation  site  of  a  crack  and  its  tip  during  propagation,  and  uses  a 
micromorphic  continuum  approach  in  the  overlap  and  adjacent  regions  to  provide  proper 
BCs  on  the  DNS  region,  as  well  as  an  overlay  continuum  to  which  to  project  the  underlying 
particle-scale  mechanical  response  (stress,  internal  state  variables  [ISVs] ) .  The 
micromorphic  continuum  constitutive  model  accounts  for  the  inherent  length-scale  of 
damaged  fracture  zone  at  the  particle-scale,  and  thus  includes  the  kinematics  to  enable  the 
proper  coupling  with  the  fractured  DNS  particle  region.  Outside  of  the  DNS  region,  a 
micromorphic  extension  of  existing  continuum  model(s),  with  the  particular  model(s)  to  be 
determined  based  on  ARL  needs,  of  material  behavior  is  used. 

This  SOW  calls  for  development  of  the  formulation  and  finite  element  implementation  of  a 
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finite  strain  micromorphic  inelastic  constitutive  model  to  bridge  particle-scale  mechanics  to 
the  continuum-scale.  The  desired  result  is  formulation  of  a  model  that  enables  more 
complete  understanding  of  the  role  of  microstructure-scale  physics  on  the 
thermomechanical  properties  and  performance  of  heterogeneous  materials  of  interest  to 
ART.  These  materials  could  include,  but  are  not  limited  to,  the  following:  ceramic 
materials,  energetic  materials,  geological  materials,  and  urban  structural  materials. 

2.1.1  Specific  Tasks 

Specific  tasks,  and  summary  of  what  was  accomplished  for  each  task. 


1.  Investigate  and  assess  specific  needs  of  ARL  researchers  with  regards  to  multiscale 
modeling  of  heterogeneous  particulate  materials.  Determine,  following  discussion  with 
ARL  materials  researchers,  the  desired  classical  continuum  constitutive  model  to  be 
reformulated  as  a  micromorphic  continuum  constitutive  model  and  used  in  the  region 
outside  and  overlapping  partially  the  DNS  window,  for  material(s)  of  interest  to 
ARL.  For  example,  polycrystalline  ceramics  models  include  those  of  Johnson  and 
Holmquist  (1999)  or  Rajendran  and  Grove  (1996)  and  energetic  materials  include 
those  following  Dienes  et  al.  (2006).  A  finite  strain  Drucker-Prager  pressure-sensitive 
elastoplasticity  model  [Regueiro,  2009]  was  selected  as  a  simple  model  approximation 
to  start,  with  future  extension  to  the  more  sophisticated  constitutive  model  forms 
mentioned  in  the  task  1.  This  model  is  presented  in  section  2.3.3. 

2.  Formulate  theory  and  numerical  algorithms  for  a  finite  strain  micromorphic  inelastic 
constitutive  model  to  bridge  particle- scale  mechanics  to  the  continuum- scale  based  on 
the  decided  constitutive  equations  from  task  1. 

See  the  summary  for  task  1. 

3.  Initiate  finite  element  implementation  of  the  formulated  finite  strain  micromorphic 
inelastic  constitutive  model  in  a  continuum  mechanics  code. 

The  finite  element  implementation  has  been  initiated  in  the  password-protected 
version  of  Tahoe  tahoe.colorado.edu,  where  the  open  source  version  is  available  at 
tahoe.cvs.  sourceforge.net.  This  report  focuses  on  the  theory;  while  details  of  the 
finite  element  implementation  and  numerical  examples  will  follow  in  journal  articles 
and  a  future  report. 

4.  Interact  with  ARL  researchers  in  order  to  improve  mutual  understanding  (i.e., 
understanding  of  both  principal  investigators  [Pis]  and  of  ARL )  with  regards  to 
dynamic  fracture  and  material  degradation  in  bound  particulate  materials  and 
associated  numerical  modeling  techniques. 


We  are  continuing  to  interact  with  ARL  researchers  regarding  their  needs  for  this 
research  problem. 


5.  Formulate  an  algorithm  to  couple  finite  strain  micromorphic  continuum  finite 
elements  to  DNS  finite  elements  of  bound  particulate  material  through  an  overlapping 
region. 

The  formulated  algorithm  is  presented  in  section  2.5. 

6.  Initiate  implementation  of  coupling  algorithm  in  task  5  using  FE  code  Tahoe  (both  for 
micromorphic  continuum  and  DNS).  Future  extension  can  be  made  for  coupling 
micromorphic  model  (Tahoe)  to  DNS  model  (ARL  or  other  FE,  or  particle/meshfree, 
code). 

The  coupling  algorithm  has  been  initiated  for  a  finite  element  and  discrete  element 
coupling.  Extension  to  other  DNS  models  of  the  grain-scale  response  is  part  of  future 
work  (see  section  2.5). 

2.2  Nonlinear  Micromorphic  Continuum  Mechanics 
2.2.1  Kinematics 

Figure  3  illustrates  the  mapping  of  the  macroelement  and  microelement  in  the  reference 
configuration  to  the  current  configuration  through  the  deformation  gradient  F  and 
microdeformation  tensor  x-  The  macroelement  continuum  point  is  denoted  by  P(X,  E)  and 
p(x,  £,  t)  in  the  reference  and  current  configurations,  respectively,  with  centroid  C  and  c. 
The  microelement  continuum  point  centroid  is  denoted  by  C'  and  d  in  the  reference  and 
current  configurations,  respectively.  The  microelement  is  denoted  by  an  assembly  of 
particles,  but  in  general  represents  a  grain/particle/fiber  microstructural  sub-volume  of  the 
heterogeneous  material.  The  relative  position  vector  of  the  microelement  centroid  with 
respect  to  the  macroelement  centroid  is  denoted  by  S  and  £(X,  S,  t)  in  the  reference  and 
current  configurations,  respectively,  such  that  the  microelement  centroid  position  vectors 
are  written  as  (figure  3)  (Eringen  and  Suhubi,  1964;  Eringen,  1999) 


X'k  —  Xk  +  ^k-,  x'k  —  Xk(X-fi)  +  £fc(X,  E,  t)  (1) 

Eringen  and  Suihubi  (1964)  assumed  that  for  sufficiently  small  lengths  ||E||  <C  1  (  ||  •  ||  is 
the  L2  norm),  (  is  linearly  related  to  E  through  the  microdeformation  tensor  y,  such  that 


£fc(X,  E,f)  —  Xfcft:(X,  t)EK 


(2) 


9 


Figure  3.  Map  from  reference  Bq  to  current  configuration  B  accounting  for  relative  posi¬ 
tion  S,  £  of  microelement  centroid  C',  c!  with  respect  to  centroid  of  macroele¬ 
ment  C,  c.  F  and  y  can  load  and  unload  independently  (although  coupled 
through  constitutive  equations  and  balance  equations),  and  thus  the  addi¬ 
tional  current  configuration  is  shown. 


where  then  the  spatial  position  vector  of  the  microelement  centroid  is  written  as 


x’k  =  xk(X,  t )  +  XkKpt,  t) ZK  (3) 

This  is  equivalent  to  assuming  an  affine,  or  homogeneous,  deformation  of  the  macroelement 
differential  volume  dV  (but  not  the  body  B ;  i.e.,  the  continuum  body  B  is  expected  to 
experience  heterogeneous  deformation  because  of  y,  even  if  BCs  are  uniform).  It  also 
simplifies  considerably  the  formulation  of  the  micromorphic  continuum  balance  equations 
as  presented  in  (Eringen,  1964;  Eringen,  1999).  This  micro  deformation  y  is  analogous  to 
the  small  strain  microdeformation  tensor  -0  in  Mindlin  (1964),  physically  described  in  his 
figure  1.  Eringen  (1968)  also  provides  a  physical  interpretation  of  y  generally,  but  then 
simplies  for  the  micropolar  case.  For  example,  y  can  be  interpreted  as  calculated  from  a 
microdisplacement  gradient  tensor  $  as  y  =  1  +  $,  where  $  is  not  actually  calculated  from 
a  micro  displacement  vector  u',  but  a  u'  can  be  calculated  once  $  is  known  (see  equation 
265).  The  microelement  spatial  velocity  vector  (holding  X  and  S  fixed)  is  then  written  as 
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v'k  =  x'k  =  Xk  +  L  =  vk  +  Vldil 


(4) 


where  ^  =  XkK^K  =  XhkXkiAi  =  uki£,h  vk  is  the  macroelement  spatial  velocity  vector, 
vki  =  XkKXx \  (u  =  XX1)  the  microgyration  tensor,  similar  in  form  to  the  velocity  gradient 

vk,i  =  FkKFxl  (i  =  FF-1). 


Now  we  take  the  partial  spatial  derivative  of  (equation  3)  with  respect  to  the  reference 
microelement  position  vector  X'K)  to  arrive  at  an  expression  for  the  microelement 
deformation  gradient  F’kJ<  as  (see  appendix  A) 


F' 

rkK 


v  mt  ^  ,  dXkL(X,t)r_ 
*kK{-^,  t)  H - — - 1 — 'Z/ 


OX 


K 


,  /  fv  17  /V  ^  ^XfeM(X,  t)  „  \  dhA 

+  XfcA(X,  t)  -  FkA(X,  t) - 


dX, 


dX 


K 


(5) 


where  the  deformation  gradient  of  the  macroelement  is  FkK  =  dxk(X.,t)/ dXj<.  The 
microelement  deformation  gradient  FkK  maps  microelement  differential  line  segments 
dx'k  =  F'kKdX' K  and  volumes  dv'  =  J'dV' ,  where  J'  =  detF'  is  the  microelement  Jacobian 
of  deformation.  This  is  presented  for  generality  of  mapping  stresses  between  £>0  and  £>,  £>o 
and  B,  B  and  B,  but  will  not  be  used  explicitly  in  the  constitutive  equations  in  section 
2.3.3. 


2.2.2  Micromorphic  Balance  Equations  and  Clausius-Duhem  Inequality 

Using  the  spatial  integral- averaging  approach  in  Eringen  and  Suhubi  (1964),  we  can  derive 
the  balance  equations  and  Clausius-Duhem  inequality  summarized  in  equation  57.  The 
rationale  of  this  integral-averaging  approach  over  dv  and  B  in  the  current  configuration  is 
to  assume  the  classical  balance  equations  in  the  microelement  differential  volume  dv'  must 
hold  over  integrated  macroelement  differential  volume  dv ,  in  turn  integrated  over  the 
current  configuration  of  the  body  in  B.  This  approach  is  applied  repeatedly  to  derive  the 
micromorphic  balance  equations  in  equation  57. 

Balance  of  mass:  The  microelement  mass  m!  over  dv  can  be  expressed  as 


m 


PodV 


where  p'0  =  p'J',  J'  =  detF'.  Then,  the  conservation  of  microelement  mass  m'  is 


(6) 
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=  [  p'dv'  =  [  p'J'dV' 

Dt  Jdv  Dt  JdV 


Thus,  the  pointwise  (localized)  balance  of  mass  over  dv  is 


Dp'  .dv', 

~M+l,Mr 


Now,  consider  the  balance  of  mass  of  solid  over  the  whole  body  B.  We  start  with  the 
integral-average  definition  of  mass  density: 


pdv  =  /  p'dv' 


The  total  mass  m  of  body  B  is  expressed  as 


rn  —  I  pdv  =  /  p  dv  =  /  p  J  dV‘ 

JB  JB  IJdv  \  Jb0  IJdV 


Then,  for  conservation  of  mass  over  the  body  B,  we  have 


B0  UdV 


Bmdv> 

Dt 


f  f  ^P'  /  dv[  , 

L  L  « +  " 


Then,  the  balance  of  mass  in  B  leads  to  the  standard  result 
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Dm 

~Dt 


D 

~Dt 


/  pdv  =  0 

Jb 

D(pJ) 


1 B0 


Dt 


-dV 


DP  ,  dvi  . 
-+p—)dv  =  0 


(12) 


Localizing  the  integral,  we  have  the  pointwise  satisfaction  of  balance  of  mass  for  a  single 
constituent  (in  this  case,  solid)  material: 


Dp  dvi 
Dt  +PdXl 


0 


(13) 


Balance  of  microinertia: 


Given  that  Zk  is  the  position  of  microelement  dV'  centroid  C'  in  the  reference 
configuration  with  respect  to  the  mass  center  of  the  macroelement  dV  centroid  C  (see 
figure  3),  we  have  the  result 


[  p'0ZKdV'  =  0  (14) 

JdV 

This  can  be  thought  of  as  the  first  mass  moment  being  zero  because  of  the  definition  ZK  as 
the  “relative”  position  of  C'  with  respect  to  C  (the  mass  center  of  dV)  (Eringen,  1999). 
The  second  mass  moment  is  not  zero,  and  in  the  process  a  microinertia  Ikl  hr  is  defined 
as 


PoIKLdV=  [  p0ZKZLdV  (15) 

JdV 

Likewise,  a  microinertia  in  B  is  dehned  as 
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pikidv 


(16) 


def 


'  dv 


p'ikiidv' 

f> XkK^KXlL^Ldv' 


'  dv 


XkKXlL  /  Po^K^L 
J  dv 


dV' 


XkKXlL  Pod KLdV  —  XkKXlL  pd K  Ldv 
=►  IKL  =  XKkXlliki 


The  balance  of  microinertia  in  Bq  is  then  defined  as 


(17) 


D 

Dt 


PoIKLdV 


/  Po  —  dV  =  0 
JBo 

DIkl 


Dt 


Dt 


Xk\Xli 


Di 


kl 


Dt 


^ka^al  Vialak 


PX.KkX.Ll 


Di 


kl 


Dt 


^ka^al  ^la^ak  )  dv  0 


(18) 


Localizing  the  integral,  and  factoring  out  PX.k\,X.li  ■>  the  pointwise  balance  of  microinertia  in 
B  is 


Pin 

~Dt 


V kcPal  V iDak  0 


(19) 


Balance  of  linear  momentum  and  the  first  moment  of  momentum:  To  derive  the 
micromorphic  balance  of  linear  momentum  and  the  first  moment  of  momentum  (different 
than  angular  momentum),  Eringen  and  Suhubi  (1964)  followed  a  weighted  residual 
approach,  where  the  point  of  departure  is  that  balance  of  linear  and  angular  momentum  in 
the  microelement  dv'  over  dv  are  satisfied: 


+  P'(fl  -  <)  =  0  (20) 

a'ik  =  °'ki  (21) 

where  microelement  Cauchy  stress  o'  is  symmetric  (macroelement  Cauchy  stress  o  will  be 
shown  to  be  symmetric).  Using  a  smooth  weighting  function  (f V  (to  be  defined  for  three 
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cases),  the  weighted  average  over  B  of  the  balance  of  linear  momentum  on  dv  is  expressed 
as 


Wm  +  P'ifk 


=  0 


(22) 


where  (•)/i  =  &(•)' /dx[.  Applying  the  chain  rule  (d'cr^.)  j  =  4>'iCr'lk  +  4>><Jikii  we  can  rewrite 
equation  22  as 


We  consider  three  cases  for  the  weighting  function  (j)'  leading  to  three  separate 
micromorphic  balance  equations  on  B: 

1.  (j)'  —  1,  balance  of  linear  momentum 

2.  q V  =  enmkx'm,  balance  of  angular  momentum,  where  enmk  is  the  permutation  tensor 
(Holzapfel,  2000) 

3.  (j)'  =  x'm,  balance  of  the  first  moment  of  momentum 

Substituting  these  three  choices  for  <fi'  into  equation  24,  we  can  derive  the  respective 
micromorphic  balance  equations  on  B: 

1.  (/>'  =  1,  balance  of  linear  momentum: 


The  spatial-averaged  definitions  of  unsymmetric  Cauchy  stress  atk,  body  force  /*,,  and 
acceleration  ak  are  used  to  derive  the  micromorphic  balance  of  linear  momentum: 


aikjiida 

def 

/  <7'lkn'lda' 

(26) 

J  da 

pfkdv 

def 

[  P'f'kdv' 

(27) 

J  dv 

pakdv 

def 

/  p  a'kdv' 

J  dv 

(28) 
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From  equation  25  and  equations  26-28,  there  results 


Id  B 


crikUida  +  /  p(fk  -  ak)dv  =  0 


+  p(fk  -  a*)]  dv  =  0 


(29) 

(30) 


Localizing  the  integral,  we  have  the  pointwise  expression  for  micromorphic  balance  of 
linear  momentum 


<?ik,i  +  p{fk  ~  ak)  =  0  (31) 

Note  that  the  macroscopic  Cauchy  stress  aik  is  unsymmetric. 

2.  <j)'  =  enmkx'm ,  x'm  =  xm  +  Cmi  balance  of  angular  momentum: 


'  8B 


/  enmkix^a'Mda'  >+  <  enmk  \-xmlo'lk  +  p'x'm(f'k  -  a'k)]  dv '  \  =  0 

Ida  J  JB  v  J dv 


d  B  yJ  da 


€"rimk  .  O mk  P  £>m)(fk  ^k) .  (  0 


B  \  J  dv 


(32) 


where  x'm  l  =  dx'm/dx[  =  5mi.  We  analyze  the  terms  in  equation  32,  using 
a'k  =  ak  +  (k  and  £k  =  (vkc  +  vkbVbc)£, c,  such  that 


^timfc((^tn  "F  £m)crlk)nlda 


I  £"nmk%m  / 
IdB 


'  da 


def  , 

=  ( Tn-Tiida 


“I”  /  ^ nmk  I 


IdB 


’  da 


def  , 

—  m,lkm^'l  ®a 


& nmk  I  \pC"m  ^Ik^l 
JOB 


& nmk  I  [^ra/c  H“  XinO' Ik. I  1Xlikm,l\  dv 

JB 


(33) 
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I  & nmk  \p  (*^m  “1“  £>m)  f k\  dv  r  I  &nmk%m  I  P  f 
'dv  J  JJ3  Jdv 

def  P  , 

=  PJkdv 


“1“  /  €"rimk  I  P  f  k^m^V 


'  B  Jdv 


d=p(-kmdv 


6 nmk  I  (jJniPfk  T  P^km)  dv 
JB 


(35) 


^ nmk  \p  (%m  “1“  ^m)(  ^k) .  d^  f 


B  KJdv 


& nmk  I  \  I  P  “1“  ^m^k  “1“  ^m^k 

’B  l  Jdv 


~^~^m^,k^)dv 


&nmk 


XmCLk  I  p'dv' 

dv 


def  , 

=  pdv 


H ~%m(Pkc  “t-  ^kb^bc)  I  P  £>cdv  ~\~CLk  I  P  ^mdv 

’  dv  J  dv 


=0 


=0 


P  £,k£,mdv 


’  dv 


def  , 

—  P^kmdv 


&nmk  I  [j’rnPJk  "1“  pjJkm]  dv 
JB 


(36) 


where  mikm  is  the  higher  order  (couple)  stress,  smk  is  the  symmetric  microstress,  £km 
is  the  body  force  couple,  and  c Jkm  is  the  micro-spin  inertia.  Combining  the  terms,  we 
have 


&nmk 


%m(&lk,l  p(fk  QJk)')  &mk  $mk  “1“  TTllkm,l  P^km  ^ km ) 


=0 


dv  =  0 


€"nmk  I  [^m/c  $mk  lkm,l  p(^km  dv  0 

Jb 


(37) 


Thus,  upon  localizing  the  integral, 
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(38) 

(39) 


Cnm.k  \Prnk  Smk  T  Wilkin, l  "I-  P^km  ^4cm)]  0 

W[mk]  S[mk\  d-‘Wll[km],l  “I"  P{^[km]  ^[kin])  0 

=0 


resulting  in 


W[mk]  4~  ^lZ[fcm],Z  T  P(^[km]  km] )  0  (40) 

where  the  antisymmetric  definition  cr[mfc]  =  ( amk  —  akm)/2.  Equation  40  is  the 
pointwise  balance  of  angular  momentum  on  B ,  providing  three  equations  to  solve  for 
a  microrotation  vector  ipk  (Eringen,  1968).  But  we  want  to  solve  for  the  general 
nine- dimensional  microdeformation  tensor  XkK,  thus  we  need  six  more  equations. 

The  balance  of  the  first  moment  of  momentum  provides  these  additional  equations. 

3.  (f)'  =  x'm,  balance  of  the  first  moment  of  momentum:  The  analysis  follows  that  used 
for  the  balance  of  angular  momentum,  except  we  do  not  multiply  by  the  permutation 
tensor  enmk .  Thus,  we  may  write  equation  38  directly  without  enmk : 

Wmk  Smk  T  Tfl lkm,l  T  P^km  ^ km )  0  (41) 

This,  in  general,  provides  nine  equations  to  solve  for  a  microdisplacement  gradient 
tensor  &kK  through  the  definition  XkK  =  $kK  +  ®kK-  We  note  that  euation  41 
encompasses  equation  40  (the  three  antisymmetric  equations),  and  provides  six 
additional  equations  (the  symmetric  part  of  equaiton  41)  (Eringen  and  Suhubi,  1964). 


Balance  of  energy:  It  is  assumed  the  classical  balance  of  energy  equation  holds  in 
microelement  dv'  over  macroelement  dv  as 


+  l't,k  +  p'r')dv ' 


(42) 


where  e'  is  the  microinternal  energy  density  per  unit  mass,  q'k  the  micro-heat  flux,  and  r' 
the  micro-heat  source  density  per  unit  mass.  This  is  then  integrated  to  hold  over  the  whole 
body  B  as 


+  qi  k  +  p'r')dv' 


(43) 
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The  individual  terms  in  equation  43  can  be  analyzed,  using  v[  =  Vi  +  ^  =  Vi  +  zym£m, 
a[  =  ai  +  6,  and  a'klk  =  p\a[  -  //): 


°'kiKkdv' 


D 

Dt 


(p0edV)  =  p0ed,V  =  pedu 


(44) 


(45) 


/  /  /  7  / 

(j  kiv^da 


'  da 


’  dv 


v'ki,kvidv' 


a 


kl 


( Vl  + 


)n'kda! 


p’(a'i  -  fl)(vi  +  vim£,m)dv' 


'  da 


'  dv 


=  Vl 


aklnkda  I  c t 


kl^rn'nkda 


-Vl 


'  da 


'da 


p'ci^dv'  +vi 


'  dv 


def  ,  def  .  def  .  def  „  . 

=  (7klnkda  =mkimnkda  =  paidv  =  pjidv 


I  P  ,mdv  Vim  I  p  ^l^mdv  I  P  fl^mdv 


J  dv 


'  dv 


'  dv 


=0 


def  , 

—  pt-^lmdv 


d=  phmdv 


q'kn'kda!  =f  qknkda 


'da 


def 


prdv 


(46) 

(47) 


Substituting  these  terms  back  into  equation  43,  we  have 


pedv  =  /  (■ viakink  +  vlmmkimnk)da  -  /  vip(ai  -  fi)dv  -  /  uimP^im  ~  l im)dv 

’  JdB  JB  JB 


+  /  qknkda  +  /  prdv 

JdB  Jb 


vJcTki  k  T  p(fl  ®z))  T  Virn(^mklm^k  T  p(f)m  ^Zm)) 


$ml  &ml 


-'rVi^Oki  +  Bim^rrikim  +  qkjk  +  pr]  dv 


Localizing  the  integral,  the  pointwise  balance  of  energy  over  B  becomes 


(48) 


PV  —  VimiSml  &ml)  T  Vi)k(Jki  +  k/irn^k‘lfYlkirn  T  qk,k  T  PV 


(49) 
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Second  Law  of  Thermodynamics  and  Clausius- Duhem  Inequality:  We  assume  the  second 
law  is  valid  in  microelement  dv'  over  dv  such  that 


fdv  P'ri'dv'd=pridv 


P^dv>  >  0 


'  dv 


d=^dv 


(50) 


Note  that  no  microtemperature  6'  is  currently  introduced  (Eringen.  1999).  Integrating  over 
B ,  localizing  the  integral,  and  multiplying  by  macrotemperature  9 ,  we  arrive  at  the 
pointwise  form  of  the  second  law  as 


pdf)  -  qKk  +  -qk9,k  -  pr  >  0 


dv  — 


^dv  >  0 
0 


(51) 

(52) 


We  derive  the  micromorphic  Clausius- Duhem  inequality  by  introducing  the  Helmholtz  free 
energy  function  ip,  and  using  the  balance  of  energy  in  eqnaiton  49.  Recall  the  definition  of 
ip  (Holzapfel,  2000),  and  its  material  time  derivative  leading  to  an  expression  for  p9f\  in 
equation  52  as 


ip  =  e  —  Or)  (53) 

ip  —  e  —  Or)  —  Of)  (54) 

pOf)  =  pe  —  pOp  —  pip  (55) 

Upon  substitution  into  equation  52  and  using  equation  49,  we  arrive  at  the  micromorphic 
Clausius-Duhem  inequality: 

-p(ip  +  1)0)  +  (Jkliypk  -  m)  +  S kiVik  +  mkimvim,k  +  \qk0,k  >  0  (56) 
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Summary  of  balance  equations:  The  equations  are  now  summarized  over  the  current 
configuration  B  as 


balance  of  mass  :  +  pvkik  —  0 

>  def  r 


Pdv  =  fdv  P'dv' 


balance  of  microinertia  :  —  vkmimi  —  vimimk  =  0 

def 


pikidv  =  fdv  p'tktidv' 


balance  of  linear  momentum  :  +  p(fk  —  ak)  =  0 

crik^da  =  fda  o-jk'n'da' 
pfkdv  ==  fdv  p’fldv’ 
pakdv  =  fdv  p'a'kdv' 


balance  of  first  moment  of  momentum  :  ami  —  smi  +  rrikim,k  +  p(£im  —  w im )  =  0  /  (57) 

7  def  p  iii 


Smldv  =  fdv  <j'mldv' 
mkimnkda  =  fda  o'kl£mn'kda! 
phmdv  =  Jdv  p'flimdv' 
pOJimdv  =  fdv  pf'iiimdv' 


balance  of  energy  :  pe  =  (ski  -  crki)isik  +  <JkiVhk 

T mklm^lm,k  T  Qk,k  T  P^ 

Clausius  -  Duhem  inequality  :  -p(ip  +  pd)  +  crki(vpk  -  uik)  +  skiuik 

T TYlklm,L/lm,k  T  ~gdkP,k  ^  0 


where  D(m)/Dt  is  the  material  time  derivative,  iki  is  the  symmetric  microinertia  tensor,  aik 
is  the  unsymmetric  Cauchy  stress,  fk  the  body  force  vector  per  unit  mass,  //  is  the  body 
force  vector  per  unit  mass  over  the  microelement,  ak  is  the  acceleration,  smi  is  the 
symmetric  microstress,  mkim  the  higher  order  couple  stress,  £im  the  body  force  couple  per 
unit  mass,  ujim  the  microspin  inertia  per  unit  mass,  e  is  the  internal  energy  per  unit  mass, 
uik  is  the  microgyration  tensor,  vpk  is  the  velocity  gradient,  vim,k  is  the  spatial  derivative  of 
the  microgyration  tensor,  qk  is  the  heat  flux  vector,  r  is  the  heat  supply  per  unit  mass,  ^  is 
the  Helmholtz  free  energy  per  unit  mass,  rj  is  the  entropy  per  unit  mass,  and  6  is  the 
absolute  temperature.  Note  that  the  balance  of  first  moment  of  momentum  is  more  general 
than  the  balance  of  angular  momentum  (or  “moment  of  momentum”  [Eringen,  1962]),  such 
that  its  skew-symmetric  part  is  the  angular  momentum  balance  of  a  micropolar  continuum 
(see  equation  40)).  Recall  that  the  Cauchy  stress  a'ml  over  the  microelement  is  symmetric 
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Figure  4.  Differential  area  of  microelement  da'  within  macroelement  da  in  current  con¬ 
figuration  B. 

because  the  balance  of  angular  momentum  is  satisfied  over  the  microelement  (Eringen  and 
Suhubi,  1964). 

Physically,  the  microstress  s  defined  in  equation  574  as  the  volume  average  of  the  Cauchy 
stress  o'  over  the  microelement,  can  be  interpreted  in  the  context  of  its  difference  with  the 
unsymmetric  Cauchy  stress  as  s  —  o  (Mindlin  (1964)  called  this  the  “relative  stress”).  This 
is  the  energy  conjugate  driving  stress  for  the  microdeformation  y  through  its  microgyration 
tensor  v  =  yy^1  in  equation  575,  and  also  the  reduced  dissipation  inequality  in  the 
intermediate  configuration  equations  95  and  98  as  S  —  S  (the  analogous  stress  difference  in 
B).  In  fact,  we  do  not  solve  for  s  or  S  directly,  but  constitutively  we  solve  for  the 
difference  s  —  o  or  S  —  S  (see  equation  118).  The  higher-order  stress  m  is  analogous  to  the 
double  stress  /i  in  Mindlin  (1964)  with  physical  components  of  microstretch,  microshear, 
and  microrotation  shown  in  his  figure  2.  For  example,  mu 2  is  the  higher  order  shear  stress 
in  the  X2  direction  based  on  a  stretch  in  the  X\  direction.  Using  the  area  average  definition 
for  mum,  we  have  m^ni  =f  (1  /da)  fda  o'^^n^da' ,  where  o'u  is  the  normal  microelement 
stress  in  the  X\  direction,  and  £2  is  the  shear  couple  in  the  X2  direction. 

2.3  Finite  Strain  Micromorphic  Elastoplasticity 

This  section  proposes  a  phenomenological  bridging-scale  constitutive  modeling  framework 
in  the  context  of  finite  strain  micromorphic  elastoplasticity  based  on  a  multiplicative 
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decomposition  of  the  deformation  gradient  F  and  microdeformation  tensor  x  into  clastic 
and  plastic  parts.  In  addition  to  the  three  translational  displacement  vector  u  degrees  of 
freedom  (DOFs),  there  are  nine  DOFs  associated  with  the  unsymmetric  microdeformation 
tensor  y  (microrotation,  microstretch,  and  microshear).  We  leave  the  formulation  general 
in  terms  of  y,  which  can  be  further  simplified  depending  on  the  material  and  associated 
constitutive  assumptions  (see  Forest  and  Sievert  (2003,  2006)).  The  Clausius-Duhem 
inequality  formulated  in  the  intermediate  configuration  yields  the  mathematical  form  of 
three  levels  of  plastic  evolutions  equations  in  either  (1)  Mandel-stress  form  (Mandcl,  1974), 
or  (2)  an  alternate  ‘metric’  form.  For  demonstration  of  the  micromorphic  elastoplasticity 
modeling  framework,  J2  flow  plasticity  and  linear  isotropic  elasticity  are  initially  assumed, 
extended  to  a  pressure-sensitive  Drucker-Prager  plasticity  model,  and  then  mapped  to  the 
current  configuration  for  semi-implicit  numerical  time  integration. 

The  formulation  presented  here  differs  from  other  works  on  finite  strain  micromorphic 
elastoplasticity  that  consider  a  multiplicative  decomposition  into  elastic  and  plastic  parts 
(Sansour  1998,  Forest  and  Sievert,  2003,  2006)  and  those  that  do  not  (Lee  and  Chen,  2003; 
Vernerey  et  al.,  2007). 

Sansour  (1998)  considered  a  finite  strain  Cosserat  and  micromorphic  plastic  continuum, 
redefining  the  micromorphic  strain  measures  (see  equation  B-)  in  appendix  B)  to  be 
invariant  with  respect  to  rigid  rotations  only,  not  also  translations.  Sansour  did  not  extend 
his  formulation  to  include  details  on  a  finite  strain  micromorphic  elastoplasticity 
constitutive  model  formulation,  as  this  report  does.  Sansour  proposed  to  arrive  at  the 
higher-order  macro-continuum  by  integrally- averaging  micro-continuum  plasticity  behavior 
using  computation.  Such  an  approach  is  similar  to  computational  homogenization,  as 
proposed  by  Forest  and  Sievert  (2006)  to  estimate  material  parameters  for  generalized 
continuum  plasticity  models.  On  a  side  note,  one  advantage  to  the  micromorphic 
continuum  approach  by  Eringen  and  Suhubi  (1964)  is  that  the  integral-averaging  of  certain 
stresses,  body  forces,  and  microinertia  terms  are  already  part  of  the  formulation.  This 
becomes  especially  useful  when  computationally  homogenizing  the  underlying 
microstructural  mechanical  response  (e.g.,  provided  by  a  microstructural  FE  or  DE 
simulation)  in  regions  of  interest,  such  as  overlapping  between  micromorphic  continuum 
and  grain/particle/fiber  representations  for  a  concurrent  multiscale  modeling  approach 
(figure  12). 

Forest  and  Sievert  (2003,  2006)  established  a  hierarchy  of  elastoplastic  models  for 
generalized  continua,  including  Cosserat,  higher  grade,  and  micromorphic  at  small  and 
finite  strain.  Specifically  with  regard  to  micromorphic  finite  strain  theory,  Forest  and 
Sievert  (2003)  follows  the  approach  of  Germain  (1973),  which  leads  to  different  stress 
power  terms  in  the  balance  of  energy  and,  in  turn,  Clausius-Duhem  inequality  than 
presented  by  Eringen  (1999).  Also,  the  invariant  clastic  deformation  measures  do  not 
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match  the  sets  89  and  B-l  proposed  by  Eringen  (1999).  Upon  analyzing  the  change  in 
square  of  microelement  arc-lengths  ( ds ')2  —  ( dS ')2  between  current  B  and  intermediate 
configurations  B  (cf.  appendix  C),  then  either  set  89  or  B-l  is  unique.  Forest  and  Sievert 
(2003,  2006)  proposed  to  use  a  mix  of  the  two  sets,  i.e.  equations  8)i,  B-l2,  and  B-l3,  in 
their  Helmholtz  free  energy  function.  When  analyzing  {ds')2  —  ( dS ')2,  they  would  also  need 
(B-l)i  as  a  fourth  elastic  deformation  measure.  As  Eringen  proposed,  however,  it  is  more 
straightforward  to  use  either  set  89  or  B-l  when  representing  elastic  deformation.  The 
report  presents  both  sets,  but  we  use  equation  89.  Mandel  stress  tensors  are  identified  in 
Forest  and  Sievert  (2003,  2006)  to  use  in  the  plastic  evolution  equations.  This  report 
presents  additional  Mandel  stresses  and  considers  also  an  alternate  “metric” -form  often 
used  in  finite  deformation  elastoplasticity  modeling. 

Vernerey  et  al.  (2007)  treated  micromorphic  plasticity  modeling  similar  to  Germain  (19730 
and  Mindlin  (1964),  which  leads  to  different  stress  power  terms  and  balance  equations  than 
in  Eringen  (1999).  The  resulting  plasticity  model  form  is  thus  similar  to  Forest  and  Sievert 
(2003),  although  does  not  use  a  multiplicative  decomposition  and  thus  does  not  assume  the 
existence  of  an  intermediate  configuration.  An  extension  presented  by  Vernerey  et  al. 

(2007)  is  to  consider  multiple  scale  micromorphic  kinematics,  stresses,  and  balance 
equations,  where  the  number  of  scales  is  a  choice  made  by  the  constitutive  modeler.  A 
multiple  scale  averaging  procedure  is  introduced  to  determine  material  parameters  at  the 
higher  scales  based  on  lower  scale  response. 

In  general,  in  terms  of  a  multiplicative  decomposition  of  the  deformation  gradient  and 
microdeformation,  as  compared  to  recent  formulations  of  finite  strain  micromorphic 
elastoplasticity  reported  in  the  literature  (just  reviewed  in  preceding  paragraphs),  we  view 
our  approach  to  be  more  in  line  with  the  original  concept  and  formulation  presented  by 
Eringen  and  Suhubi  (1964),  Eringen  (1999),  which  provide  a  clear  link  between 
microelement  and  macroelement  deformation,  balance  equations,  and  stresses.  Thus,  we 
believe  our  formulation  and  resulting  elastoplasticity  model  framework  is  more  general 
than  what  has  been  presented  previously.  The  paper  by  Lee  and  Chen  (2003)  also  follows 
closely  Eringen’s  micromorphic  kinematics  and  balance  laws,  but  does  not  treat 
multiplicative  decomposition  kinematics  and  subsequent  constitutive  model  form  in  the 
intermediate  configuration,  as  this  report  does.  We  demonstrate  the  formulation  for  three 
levels  of  J-2  plasticity  and  linear  isotropic  elasticity,  as  well  as  pressure-sensitive 
Drucker-Prager  plasticity,  and  numerical  time  integration  by  a  semi-implicit  scheme  in  the 
current  configuration  B. 

2.3.1  Kinematics 

We  assume  a  multiplicative  decomposition  of  the  deformation  gradient  (Lee,  1969)  and 
microdeformation  (Sansour,  1998;  Forest  and  Sievert,  2003,  2006)  (figures  5),  such  that 
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F  _  pepP 


X  =  xexp 


(58) 


FkK  ^kK^KK  ’  XkK  XkxXxK 

Given  the  multiplicative  decompositions  of  F  and  y,  the  velocity  gradient  and 
microgyration  tensors  can  be  expressed  as 


t 

Vl,k 

V 

Vlk 


FF  *  —  pcpc 


—  ■ 


he  ne-1  ,  ne  jp  ne- 1  ne  , 

^ IA r  Ak  '  Ck  *~lk  + 

Lp-  -  =  F-  Fp_1- 

BC  BB  BC 


'Ik 


XX  1  =  *V  1  +  XeLx,PX 


,e  —  1 


=  ue  +  ul 

,p 


XaX'Ak  +  XlBl&PcXecl  =  vik  +  "lk 


LX-’P  =  y  P  YP 
BC  '-DD-V 


(59) 


(60) 


In  the  next  section,  the  Clausius-Duhem  inequality  requires  the  spatial  derivative  of  the 
microgyration  tensor,  which  split  into  elastic  and  plastic  parts  based  on  (60).  Thus,  it  is 
written  as 


Vv 

=  Vz/e  +  Vz/p 

(61) 

L'lrrijk 

^Im^k  ^lm,k 

V fm,k  =  XlA,kXe~Am  ~  vtnXnDjXt™ 

(62) 

Ulm,k  —  (xic,k  Xqa  XlE  X^A,k  ~  X-IF  ^  FG  X-GA,k)  XAtu 

~p'laXaA,kX  Am 

(63) 

The  spatial  derivative  of  the  elastic  microdeformation  tensor  Vye  is  analogous  to  the  small 
strain  microdeformation  gradient  K  in  Mindlin  (1964),  and  its  physical  interpretation  in 
figure  2  of  Mindlin  (1964).  For  example,  X112  is  an  elastic  microshear  gradient  in  the  X2 
direction  based  on  a  microstretch  in  the  X\  direction.  Furthermore,  just  as  differential 
macroelement  volumes  map  as 


dv  =  JdV  =  JeJpdV  =  JedV 


(64) 


where  Je  =  det  Fe  and  Jp  =  det  Fp,  then  microelement  differential  volumes  map  as 
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Figure  5.  Multiplicative  decomposition  of  deformation  gradient  F  and  microdeformation 
tensor  x  into  clastic  and  plastic  parts,  and  the  existence  of  an  intermediate 
configuration  B.  Since  Fe.  Fp,  y6,  and  xp  can  load  and  unload  independently 
(although  coupled  through  constitutive  equations  and  balance  equations),  ad¬ 
ditional  configurations  are  shown.  The  constitutive  equations  and  balance 
equations  presented  in  the  report  govern  these  deformation  processes,  and  so 
generality  is  preserved. 
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dv'  =  j'dv'  =  je,jp,dv'  =  r'dV 


(65) 


where  Je/  =  detFe/  and  Jpl  =  detFp/.  Fe/  and  Fpl  have  not  been  defined  from  equation  5), 
and  are  not  required  for  formulating  the  final  constitutive  equations.  Likewise,  according  to 
micro-  and  macroelement  mass  conservation,  mass  densities  map  as 


p0  =  pJ  =  pJeJp  =  pjp  (66) 

p’Q  =  p'j'  =  p'r'jv  =  p'jp'  (67) 


This  last  result  was  achieved  by  using  a  volume- average  definition  relating  macroelement 
mass  density  to  microelement  mass  density  as 


pdv  '=  f  p' dv'  ,  PodV  =f  f  p'0dV'  ,  pdV  =f  f  p'dV'  (68) 

Jdv  JdV  JdV 

This  volume  averaging  approach  by  Eringen  and  Suhubi  (1964)  is  used  extensively  in 
formulating  the  balance  equations  and  Clausius-Duhem  inequality  in  section  2.2.2. 

2.4  Clausius-Duhem  Inequality  in  B 

This  section  focusses  on  the  Clausius-Duhem  inequality  mapped  to  the  intermediate 
configuration  to  identify  evolution  equations  for  various  plastic  deformation  rates  that 
must  be  defined  constitutively,  and  their  appropriate  conjugate  stress  arguments  in  B. 

From  a  materials  modeling  perspective,  it  is  often  preferred  to  write  the  Clausius-Duhem 
inequality  in  the  intermediate  configuration  B,  which  is  considered  naturally  elastically 
unloaded,  and  formulate  constitutive  equations  there.  The  physical  motivation  lies  with 
earlier  work  by  Kondo  (1952),  Bilby  et  al.  (1955),  Kroner  (1960),  and  others,  who  viewed 
dislocations  in  crystals  as  defects  with  associated  local  clastic  deformation,  where 
macroscopic  elastic  deformation  could  be  applied  and  removed  without  disrupting  the 
dislocation  structure  of  a  crystal.  More  recent  models  extend  this  concept,  such  as  papers 
by  Clayton  et  al.  (2005,  2006)  and  references  therein.  The  intermediate  configuration  B 
can  be  considered  a  “reference”  material  configuration  in  which  fabric/texture  anisotropy 
and  other  inelastic  material  properties  can  be  defined.  Thus,  details  on  the  mapping  to  B 
are  given  in  this  section.  Recall  that  the  Clausius-Duhem  inequality  in  equation  57e  was 
written  using  localization  of  an  integral  over  the  current  configuration  B,  such  that 
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(69) 


p{fp  +  V6)  +  &ki(vi,k  ~ 


V Ik )  “1“  $kl^lk  F  ^T'klm^/lm,k  “1“  gQk@,k 


dv  >  0 


Using  the  microelement  Piola  transform  a'ki  =  Fel kRS'RkFel ti/  Je'  and  Nanson’s  formula 
n'kda!  =  Jel  Fj?-1  N'^dA' ,  the  following  mappings  of  the  area-averaged  unsymmetric  Cauchy 
stress  er,  volume-averaged  symmetric  microstress  s,  and  area-averaged  higher  order  couple 
stress  m  terms  are  obtained  as 
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—  mklm 


where  S'^L  is  the  symmetric  second  Piola-Kirchhoff  stress  in  the  microelement  intermediate 
configuration  (over  dV),  Spp  is  the  unsymmetric  second  Piola-Kirchhoff  stress  in  the 
intermediate  configuration  B ,  S ppp  is  the  symmetric  second  Piola-Kirchhoff  microstress  in 
the  intermediate  configuration  B ,  Mklm  is  the  higher  order  couple  stress  written  in  the 
intermediate  configuration,  and  Nfp  the  unit  normal  on  dA.  In  general,  Fe/  ^  Fe,  but  the 
constitutive  equations  in  section  2.3.3  do  not  require  that  Fe/  be  defined  or  solved. 

Using  the  mappings  for  p  and  dv,  and  the  Piola  transform  on  q the  Clausius- Duhem 
inequality  can  be  rewritten  in  the  intermediate  configuration  as 
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Individual  stress  power  terms  in  equation  73  can  be  additively  decomposed  into  clastic  and 
plastic  parts  based  on  equations  59-61.  Using  equation  61,  the  higher  order  couple  stress 
power  can  be  written  as 
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where  the  spatial  derivative  with  respect  to  the  intermediate  configuration  B  can  be 
defined  as 


def 
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The  other  stress  power  terms  using  equations  59  and  60  are  written  as 
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where  =  Fe^Fe^  is  the  right  elastic  Cauchy-Green  tensor  Ce  =  FeI  Fe  in  B ,  and 
=  FeiixeiE  an  elastic  deformation  measure  in  B  as  \j/e  =  FeTye  (cf.  appendix  C). 

Similar  to  Eringen  and  Suhubi  (1964)  for  a  micromorphic  clastic  material,  the  Helmholtz 
free  energy  function  in  B  is  assumed  to  take  the  following  functional  form  for  micromorphic 
elastoplasticity  as 


P^(Fe,  Xe,  Vye,  Z,  Zx,  VZX,  6)  (79) 

Xlro  xta,K .  Zk,  ZJ,  Zl L,  0) 

where  Zf>  is  a  vector  of  macro  strain-like  ISVs  in  B,  ZX  is  a  vector  of  micro  strain- like 
ISVs,  and  ZX  -f  is  a  spatial  derivative  of  a  vector  of  micro  strain-like  ISVs.  Then,  by  the 
chain  rule 
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;‘free  energy  per  unit  mass”  assumption  is  that 


D('j^i  )  =  p^  +  p^  =  ~(p^)j^  +  P'lp  =>  P'lp  =  (P^)j-  +  — (81) 

where  we  used  the  result  p  =  D(p0/ Jp)/ Dt  =  —pJp/Jp.  Substituting  equations  74-78  and 
equations  80,  81  into  equation  73,  and  using  the  Coleman  and  Noll  (1963)  argument  for 
independent  rate  processes  (independent  xt^,  D(xekM,R)/ Dt,  and  6),  the 
Clausius-Duhem  inequality  is  satisfied  if  the  following  constitutive  equations  hold: 
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(84) 

(85) 


For  comparison  to  the  result  reported  in  equation  6.3  of  Eringen  and  Suhubi  (1964),  we 
map  these  stresses  to  the  current  configuration,  using 
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The  equations  match  those  in  (6.3)  of  Eringen  and  Suhubi  (1964)  if  elastic,  i.e.  Fe  =  F. 
Xe  =  x-  We  prefer,  however,  to  express  the  Helmholtz  free  energy  function  in  terms  of 
invariant — with  respect  to  rigid  body  motion  on  the  current  configuration  B — elastic 
deformation  measures,  such  as  the  set  proposed  by  Eringen  and  Suhubi  (1964)  as 


C'KL  —  ^kK-^kL  >  ^KL  —  Fj kh'XkL  >  ^ KLM  ~  ^ kKXkL,M  (^9) 

We  have  good  physical  interpretation  of  Fe  (and  Fp)  from  crystal  lattice  mechanics  (Bilby 
et  al.,  1955;  Kroner,  1960;  Lee  and  Liu,  1967;  Lee,  1969),  while  the  elastic  microdeformation 
Xe  has  its  interpretation  in  figure  5  of  this  report  (elastic  deformation  of  microelement)  and 
also  figure  1  of  Mindlin  (1964)  for  small  strain  theory.  The  spatial  derivative  of  elastic 
microdeformation  Vye  has  it  physical  interpretation  in  figure  2  of  Mindlin  (1964),  and  was 
earlier  in  this  report  described  for  example,  as  xh  2  is  the  niicroshear  gradient  in  the  X2 
direction  based  on  a  stretch  in  the  x\  direction  (although  directions  are  not  exact  here 
because  of  the  spatial  derivative  with  respect  to  the  intermediate  configuration  &).  The 
Helmholtz  free  energy  function  ^  per  unit  mass  is  then  written  as 


Ce,^e,Te,Z,Zx,VZx,9)  (90) 

pWii,  n-i.  n tiM,  zr,  %  %,  o) 

and  the  constitutive  equations  for  stress  result  from  equations  82-84  as 
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where  sym  [•]  denotes  the  symmetric  part.  These  stress  equations  91-93  when  mapped  to 
the  current  configuration  are  the  same  as  equations  6.9-11  in  Eringen  and  Suhubi  (1964)  if 
there  is  no  plasticity,  i.e.  Fe  =  F  and  \e  =  y.  To  consider  another  set  of  elastic 
deformation  measures  and  resulting  stresses,  refer  to  appendix  B. 

The  thermodynamically-conjugate  stress-like  ISVs  are  defined  as 


O/  —  /  \ 


(94) 


which  are  used  in  the  evolution  equations  for  plastic  deformation  rates,  as  well  as  multiple 
scale  yield  functions,  where  we  assume  scalar  Z,  Zx,  VZX,  and  Q,  Qx ,  QVx.  The  stress-like 
ISVs  in  section  2.3.3  are  physically  interpreted  as  yield  stress  Q  and  Qx  for 
macro-plasticity  (stress  S  calculated  from  elastic  deformation)  and  micro-plasticity  (stress 
difference  E  —  S  calculated  from  elastic  deformation),  respectively,  while  QVx  is  a  higher 
order  yield  stress  for  microgradient  plasticity  (higher  order  stress  M  calculated  from 
gradient  elastic  deformation). 

The  remaining  terms  in  the  Clausius- Duhem  inequality  lead  to  the  reduced  dissipation 
inequality  expressed  in  localized  form  in  two  ways:  (1)  Mandel  form  with  Mandel-like 
stresses  (Mandel,  1974),  and  (2)  an  alternate  ‘metric’  form.  Each  leads  to  different  ways  of 
writing  the  plastic  evolution  equations,  and  stresses  that  are  used  in  these  evolution 
equations.  From  equation  73,  the  reduced  dissipation  inequality  in  Mandel  form  is  written 
as 
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and  the  spatial  derivative  of  the  micro-scale  plastic  velocity  gradient  is 
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The  Mandel  stresses  are  SxsCgL,  CgfF]F'xg(F^g  —  Sab)F%i,  and  Mrlm^zdi  where  the 
hrst  one  is  well  known  as  the  “Mandel  stress,”  whereas  the  second  and  third  are  the 
relative  micro- Mandel-stress  and  the  higher  order  Mandel  couple  stress,  respectively.  We 
rewrite  the  reduced  dissipation  inequality  equatin  95  in  an  alternate  “metric”  form  as 


-«')#  +  IQrO,k  -  Qi<Zr  -  QxRzx  -  <?g^(%) 
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Form  of  plastic  evolution  equations:  Based  on  equation  95,  in  order  to  satisfy  the  reduced 
dissipation  inequality,  we  can  write  plastic  evolution  equations  to  solve  for  FgR,  )(PgK,  and 
XPpcK  i  in  Mandel  stress  form  as 
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solve  for  \rRf;J  and  xllt,L  =  (XkK,l  ~  x'ti\'\K  j)s\'R 

where  the  arguments  in  parentheses  (•)  denote  the  Mandel  stress  and  stress-like  ISV  to  use 
in  the  respective  plastic  evolution  equation,  where  H.  Hx,  and  HVx  denote  tensor 
functions  for  the  evolution  equations,  chosen  to  ensure  that  convexity  is  satisfied,  and  the 
dissipation  is  positive.  This  can  be  seen  for  the  evolution  equations  in  102-104  by  the 
constitutive  definitions  in  equations  120,  124,  and  128  in  terms  of  stress  gradients  of 
potential  functions  (i.e.,  the  yield  functions  for  associative  plasticity).  In  an  alternate 
“metric”  form,  from  equation  98,  we  can  solve  for  the  plastic  deformation  variables  as 


ClBVm  =  H-lR  (S,  Q)  (102) 

solve  for  FrRK  and  F‘R  =  FkKFpRlR 
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n tFga.K  -  *hs*"  WdcFFruk]  =  »lXMK  (M.QV") 

(104) 

solve  for  Xrk  L  and  xIrrl  =  (XtK,L  -  xIaXak  ^xFk 

We  use  this  “metric”  form  in  defining  evolution  equations  in  section  2.3.3. 

Remark  1.  The  reason  that  we  propose  the  third  plastic  evolution  equation  101  or  104  to 
solve  for  XPRk  i  directly  (not  calculating  a  spatial  derivative  of  the  tensor  Xrk  fr°m  a  FE 
interpolation  of  xp )  is  to  potentially  avoid  requiring  an  additional  balance  equation  to  solve 
in  weak  form  by  a  nonlinear  FE  method  (refer  to  Regueiro  et  al.  (2007)  and  references 
therein).  With  future  FE  implementation  and  numerical  examples,  we  will  attempt  to 
determine  whether  equations  101  or  104  leads  to  an  accurate  calculation  of  XPRk  l  •  In 
Regueiro  (2010),  a  simpler  anti-plane  shear  version  of  the  model  demonstrates  the  two 
ways  for  calculating  Vyp,  either  by  an  evolution  equation  like  in  equations  101  or  104,  or  a 
FE  interpolation  for  xp  and  corresponding  gradient  calculation  Vyp.  Note  that  in  Forest 
and  Sievert  (2003),  for  their  equation  1553,  they  also  propose  a  direct  evolution  of  a 
gradient  of  plastic  microdeformation. 
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2.4.1  Constitutive  Equations 


The  constitutive  equations  for  linear  isotropic  elasticity,  J2  flow  associative  plasticity,  and 
Drucker-Prager  nonassociative  plasticity  with  scalar  ISV  hardening/softening  are 
formulated.  We  define  a  specific  form  of  the  Helmholtz  free  energy  function,  yield 
functions,  and  evolution  equations  for  ISVs,  and  then  conduct  a  semi-implicit  numerical 
time  integration  presented  in  section  2.3.4. 

2.4.2  Linear  Isotropic  Elasticity  and  J2  Flow  Isochoric  Plasticity 

Helmholtz  free  energy  and  stresses:  Assuming  linear  elasticity  and  linear  relation 
between  stress-like  and  strain-like  ISVs,  a  quadratic  form  for  the  Helmholtz  free  energy 
function  results  as 
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Note  that  the  ISVs  are  scalar  variables  in  this  model,  which  are  related  to  scalar  yield 
strength  of  the  material  at  two  scales,  macro  and  micro,  and  H  and  Hx  are  scalar 
hardening/softening  parameters,  and  is  a  symmetric  second  order  hardening/softening 
modulus  tensor,  which  we  assume  is  isotropic  as  =  (HVx)Skl-  Elastic  strains  are 
defined  as  (Suhubi  and  Eringen,  1964)  2 =  CJ^j  —  S^l  and  =  ~  TAe 

clastic  moduli  are  defined  for  isotropic  linear  elasticity,  after  manipulation  of  equations  in 
Suhubi  and  Eringen  (1964)  as 
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where  Arijrr  and  DrlmR  have  major  and  minor  symmetry,  while  Brzmr  and  Crirjrpq 
have  only  major  symmetry,  and  the  elastic  parameters  are  A,  n,  rj,  r,  k,  is,  a,  T\ . . .  Tn. 
Note  that  the  units  for  Ti . . .  Tn  are  stress x length2  (e.g.,  Pa.m2),  thus  there  is  a  built-in 
length-scale  to  these  clastic  parameters  for  the  higher  order  stress.  The  clastic  modulus 
tensors  Ari^r,  Brlmn,  and  Drlmn  are  not  the  same  as  in  (Eringen,  1999)  because 
different  elastic  strain  measures  were  used,  but  the  higher  order  clastic  modulus  tensor 
Crlmnpq  is  the  same.  Note  that  A  is  the  typical  linear  isotropic  elastic  tangent  modulus 
tensor,  and  A  and  fi  are  the  Lame  parameters.  After  some  algebra  using  equations  91-94, 
and  105,  it  can  be  shown  that  the  stress  constitutive  relations  are 


Brl 

—  AKLMNE%jR  +  DkbmnSrir 

+  {DrbMnEr[R  +  BremrErir)  \PL\S\p  +  &lb\ 

iP  fe  /Se-lpe 

k^’KBCNPQ1  NPQ^’LQ  1  QBC 
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^ RL 

—  ArImrErjr  +  DrbmrSrjr 

+2sym  {{DrlmrErjr  +  Brbmr£rir)  [^ia^Ab  +  ^ lb] 
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KlM 

—  CrLMNPqT RpQ 

(112) 

Q 

=  HZ 

(113) 

Qx 

=  HXZX 

(114) 

Ql 

=  HVxZ\ 

(115) 

Note  that  the  units  for  HVx  are  stress x length2  (e.g.,  Pa.m2),  thus  there  is  a  built-in 
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length-scale  to  this  hardening/softening  parameter  for  the  higher-order  stress-like  ISV. 
Assuming  elastic  deformations  are  small,  we  ignore  quadratic  terms  in  equations  110  and 
111  relative  to  the  linear  terms,  leading  to  the  simplified  stress  constitutive  equations  for 
SKL  and  as 


5 'KL 

—  {Aklmn  +  Dft :lmn)E^^  +  (-Bjv&Miv  +  :Lmn)^xin 

=  (A  +  T)(Ee^^)8KL  +  2(/x  +  v)EeRL 

(116) 

M  )  L  +  k^KL  +  u^LK 

^ KL 

—  (A  +  +  2(/i  +  cr)EeRL 

(117) 

2sym  +  ^ik] 

Note  that  the  stress  difference  used  in  equation  103  then  becomes 

S-S  =  nBeT  +  vBe  (118) 

^RL  —  Srl  =  +  U^RL 

Yield  functions  and  evolution  equations:  In  this  discussion,  three  levels  of  plastic 
yield  functions  are  defined  based  on  the  three  conjugate  stress-plastic-power  terms 
appearing  in  the  reduced  dissipation  inequality  equation  98,  with  the  intent  to  define  the 
plastic  deformation  evolution  equations  such  that  equation  98  is  satisfied.  This  allows 
separate  yielding  and  plastic  deformation  at  two  scales  (micro  and  macro)  including  the 
gradient  deformation  at  the  micro-scale.  If  only  one  yield  function  were  chosen  to  be  a 
function  of  all  three  stresses  (S,  S,  M),  then  yielding  at  the  three  scales  would  occur 
simultaneously,  a  representation  we  feel  is  not  as  physical  as  if  the  scales  can  yield  and 
evolve  separately  (although  coupled  through  balance  equations  and  stress  equations  for  S 
and  S).  Recall  the  plastic  power  terms  in  equation  98  come  naturally  from  the  kinematic 
assumptions  F  =  Ff:  Fp  and  y  =  XeXPi  and  from  the  Helmholtz  free  energy  function 
dependence  on  the  invariant  elastic  deformation  measures  C e,tye,Te,  and  the  plastic 
strain-like  ISVs  Z,  Zx,  and  V  Zx. 

macroscale  plasticity:  For  macroscale  plasticity,  we  write  the  yield  function  F  as 


F( S,a)  =  || devS ||  —  a  <  0  (119) 

||  devS ||  =  yj (devS)  :  (devS) 

(devS)  :  (devS)  =  (devS/j)  (devS/j) 

devS„  S  Sj]  -  (lo’ABSAB\  qj 1 
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where  a  is  the  macro  yield  strength  (i.e.,  stress-like  ISV  Q  =  a) 

The  definitions  of  the  plastic  velocity  gradient  Lp  and  strain-like  ISV  then  follow  as 
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(121) 

HZ 

(122) 

where  7  is  the  macroplastic  multiplier. 

Micro-scale  plasticity:  For  micro-scale  plasticity,  we  write  the  yield  function  Fx  as 


FY(£  -  S,  ax)  d=  ||dev(£  -  S)||  -  d*  <  0 


(123) 


dev(S/j  -  STJ)  d=f  (Sjj  -  Srj) 


3 CeAS&AB  -  Sab ) 


rie— 1 

LIJ 


where  ax  is  the  micro  yield  strength  (stress-like  ISV  Qx  ==f  ax).  Note  that  at  the 
microscale,  the  yield  strength  ax  can  be  determined  separately  from  the  macroscale 
parameter  a. 

Remark  2.  We  use  the  same  functional  forms  for  macro  and  micro  plasticity  ( Fx  with 
similar  functional  form  as  F,  but  different  IS  Vs  and  parameters),  but  this  is  only  for  the 
example  model  presented  here.  It  is  possible  for  the  functional  forms  to  be  different  when 
representing  different  phenomenology  at  the  micro  and  macroscales.  More  micromechanical 
analysis  and  experimental  data  are  necessary  to  determine  the  microplasticity  functional 
forms  in  the  future. 

The  definitions  of  the  microscale  plastic  velocity  gradient  Lx,p  and  strain-like  ISV  then 
follow  as 
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where  yx  is  the  micro  plastic  multiplier. 

Microscale  gradient  plasticity:  For  microscale  gradient  plasticity,  we  write  the  yield 
function  FVx  as 


FVx(M,  aVx) 


def 


||devM||  —  ||aVx||  <  0 
d evM„g  =  M-lJk  -  [Ce-1)jj 


Mabk 


(127) 


where  aVx  is  the  microgradient  yield  strength  (stress-like  ISV  QVx  =f  aVx)  Note  that  at 
the  gradient  microscale,  the  yield  strength  can  be  determined  separately  from  the  micro 
and  macroscale  parameters,  which  is  a  constitutive  assumption.  The  definitions  of  the 
spatial  derivative  of  microscale  plastic  velocity  gradient  VLX,P  and  strain-like  ISV  then 
follow  as 


\ue.  _  r  vp 
^  ldl  dm, k 


-  2#L>*w 
dFVx  d  evMjtifi 


e-ipe  _  1  cl!if  VVX 
CF  FMKi  ~  I 
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(128) 


9Mklm  |  devM 

D^,a)  clef  4..V1  x  aAX 

Dt  1  dajx  °  J||dVx|| 

(129) 

ajx  =  HVxZxl 

(130) 

where  yVx  is  the  microplastic  gradient  multiplier. 
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Remark  3.  The  main  advantage  to  defining  constitutively  the  evolution  of  the  spatial 
derivative  of  the  microscale  plastic  velocity  gradient  VLX,P  in  equation  128  separate  from 
the  microscale  plastic  velocity  gradient  Ld,p  in  equation  124  (i.e.,  no  PDE  in  XP^K)  is  to 
avoid  FE  solution  of  an  additional  balance  equation  in  weak  form.  One  could  allow  VLd,p 
and  VZX  to  be  defined  as  the  spatial  derivatives  of  Ld,p  and  Zx,  respectively,  but  then  the 
plastic  evolution  equations  are  PDEs  and  require  coupled  FE  implementation  (such  as  in 
Regueiro  et  al.  (2007)).  We  plan  to  implement  the  model,  after  time  integration  in  section 
2.3.4,  within  a  coupled  finite  element  formulation  for  the  coupled  balance  of  linear  and  first 
moment  of  momentum,  and  thus  avoiding  another  coupled  equation  to  include  in  the  FE 
equations  is  desired. 

Remark  4.  With  these  evolution  equations  in  B,  equation  120  can  be  integrated 
numerically  to  solve  for  Fp  and  in  turn  Fe,  equation  124  can  be  integrated  numerically  to 
solve  for  xp  and  in  turn  ye,  and  equation  128  can  be  integrated  numerically  to  solve  for 
Vyp  and  in  turn  Vye.  Then,  the  stresses  S,  £  —  S,  and  M  can  be  calculated  and  mapped 
to  the  current  configuration  to  update  the  balance  equations  for  FE  nonlinear  solution. 
Such  numerical  time  integration  is  carried  out  in  section  2.3.4;  and  FE  implementation  is 
ongoing  work. 

2.4.3  Drucker-Prager  Pressure-Sensitive  Plasticity 

Following  the  “metric”  form  in  equation  98,  the  J2  flow  plasticity  model  is  generalized  to 
include  pressure-sensitivity  of  yield  and  volumetric  plastic  deformation  (dilation  only  for 
now,  i.e.,  no  cap  on  the  yield  and  plastic  potential  surfaces  that  allows  plastic  compaction, 
e.g.,  pore  space  collapse). 

macroscale  plasticity:  For  macroscale  plasticity,  we  write  yield  F  and  plastic  potential  G 
functions  as 


F(S,  c)  =  || devS ||  -  (A^c  -  B^p)  <  0 

A*  =  (3*  cos  (f>  ,  B4*  =  /^sin  0  ,  (3* 


2^6 

3  +  (3  sin  cf) 


|| devS  ||  =  \J (devS)  :  (devS) 

(devS)  :  (devS)  =  (dev£/j)  (devSjj) 

devS/j  =  S,J  -  (1 0\BSAB]  Off 
G(S,  c)  =  ||devS||  —  [A^c—  B^p) 


(131) 


(132) 
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where  c  is  the  macro  cohesion,  <f>  is  the  macro  friction  angle,  i/j  is  the  macro  dilation  angle, 
and  —  1  <  /3  <  1  (/3  —  1  causes  the  Drucker-Prager  yield  surface  to  pass  through  the 
triaxial  extension  vertices  of  the  Mohr-Coulomb  yield  surface,  and  (3  —  —  1  is  the  triaxial 
compression  vertices).  Functional  forms  of  A ^  and  B ^  are  similar  to  A ^  and  B^, 
respectively,  except  (j)  is  replaced  with  -0.  The  yield  and  plastic  potential  functions  have  the 
usual  functional  form  for  pressure-sensitive  plasticity  with  cohesive  and  frictional  strength, 
as  well  as  dilatancy  (Desai  and  Siriwardane,  1984). 

Remark  5.  To  satisfy  the  reduced  dissipation  inequality,  it  can  be  shown  that  7  >  A 
(Vermeer  and  de  Borst,  1984)  which  also  has  been  verified  experimentally.  We  note  that 
cj)  >  if)  leads  to  nonassociative  plasticity,  which  violates  the  principle  of  maximum  plastic 
dissipation  (Lubliner  1990),  but  does  not  violate  the  reduced  dissipation  inequality.  It  is 
well  known  that  frictional  materials  like  concrete  and  rock  exhibit  nonassociative  plastic 
flow,  and  thus  such  features  are  also  included  here.  An  associative  flow  rule  is  reached 
when  the  friction  and  dilation  angles  are  equal,  (j)  —  i/j. 

The  definitions  of  the  plastic  velocity  gradient  Lp  and  strain-like  ISV  then  follow  as 
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where  7  is  the  macroplastic  multiplier,  and  the  stress-like  ISV  is  Q  A  c. 

Remark  6.  Note  that  the  functional  forms  of  the  plastic  evolution  equations  are  similar  to 
those  dictated  by  the  principle  of  maximum  plastic  dissipation  [Lubliner,  1990],  except  that 
a  plastic  potential  function  G  is  used  in  place  of  the  yield  function  F  (i.e.,  nonassociative). 
For  purposes  of  discussion,  we  show  the  evolution  equations  for  small  strain  plasticity: 


■p  def  .  dg 
6  =7^ 


l  def  .  dg 


(136) 


where  ep  is  the  plastic  strain,  and  (  the  strain-like  ISV.  nonassociative  plasticity  violates  the 
principle  of  maximum  plastic  dissipation,  but  we  use  similar  functional  forms  that  satisfies 
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the  reduced  dissipation  inequality  (i.e.,  the  second  law  of  thermodynamics  is  satisfied). 

Microscale  plasticity:  For  microscale  plasticity,  we  write  the  yield  Fx  and  plastic  potential 
Gx  functions  as 

Fx(£-S,cx) 


GX(S  —  S,  cx) 

where  cx  is  the  microcohesion,  0X  is  the  microfriction  angle,  -0X  is  the  micro  dilation  angle, 
and  —  1  <  (3X  <  1,  which  are  material  parameters  for  the  microscale.  Functional  forms  of 
Ax G  and  Bx G  are  similar  to  Ax ^  and  Bx respectively,  except  <f>x  is  replaced  with  ^x. 
Note  that  at  the  microscale,  the  cohesion,  friction,  and  dilation  angles  can  be  determined 
separately  from  the  macroscale  parameters,  and  likewise  the  yielding  and  plastic 
deformation. 

Remark  7.  We  use  the  same  functional  forms  for  macro  and  microplasticity  (Fx  and  Gx 
with  similar  functional  form  as  F  and  G),  but  this  is  only  for  the  example  model  presented 
here.  It  is  possible  for  the  functional  forms  to  be  different  when  representing  different 
phenomenology  at  the  micro  and  macroscales.  More  micromechanical  analysis  and 
experimental  data  are  necessary  to  determine  the  pressure-sensitive  microplasticity 
functional  forms  in  the  future. 

The  definitions  of  the  microscale  plastic  velocity  gradient  Lx,p  and  strain-like  ISV  then 
follow  as 
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(139) 


(140) 

(141) 


where  yx  is  the  microplastic  multiplier,  and  Qx  =  cx. 

Microscale  gradient  plasticity:  For  microscale  gradient  plasticity,  we  write  the  yield  FVx 
and  plastic  potential  GS7x  functions  as 


FVx( M,cVx)  =  ||devM||  -  (kLVx’^||cVx||  -  5Vx^||pVx||)  <  0 

A Vx’*  =  pVx'<t>  cos  0Vx  ,  BVx^  =  /3Vx^sin0Vx  ,  (3Vx^ 

d evMjjx  =  Mjjk  ~  Gjjlp^A 

VX  define  yr  -  - - 
Pr  —  3° abmabk 

GVx(M,cVx)  =  ||devM||  M(AVx^\\cVx\\  -  RVx^||pVx||)  (143) 

where  cVx  is  the  microgradient  cohesion,  0Vx  the  microgradient  friction  angle,  ^Vx  the 
microgradient  dilation  angle,  and  —  1  <  /3Vx  <  1,  which  are  material  parameters  for  the 
gradient  microscale.  Functional  forms  of  kLVx’^  and  BVx ^  are  similar  to  AVx’^  and  BVx,<^, 
respectively,  except  0Vx  is  replaced  with  x .  Note  that  at  the  gradient  microscale,  the 
cohesion,  friction,  and  dilation  angles  can  be  determined  separately  from  the  micro  and 
macroscale  parameters,  which  is  a  constitutive  assumption.  The  definitions  of  the  gradient 
of  microscale  plastic  velocity  gradient  VLX,P  and  strain-like  ISV  then  follow  as 
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where  yVx  is  the  micro  plastic  gradient  multiplier. 

Remark  8.  With  these  evolution  equations  in  B,  equation  133  can  be  integrated 
numerically  to  solve  for  Fp  and  in  turn  Fe,  equation  139  can  be  integrated  numerically  to 
solve  for  xp  and  in  turn  ye,  and  equation  144  can  be  integrated  numerically  to  solve  for 
Vyp  and  in  turn  Vye.  Then,  the  stresses  S,  £  —  S,  and  M  can  be  calculated  and  mapped 
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to  the  current  configuration  to  update  the  balance  equations  for  finite  element  nonlinear 
solution.  Such  numerical  time  integration  are  carried  out  in  section  2.3.4  for  the  form  of 
the  constitutive  equations  after  mapping  to  the  current  configuration. 

Mapping  constitutive  equations  to  current  configuration  B:  often,  the  constitutive 
equations  in  the  intermediate  configuration  are  mapped  to  the  current  configuration 
Eringen  and  Suhubi  [1964],  and  the  material  time  derivative  is  taken  to  obtain  an  objective 
stress  rate  and  corresponding  stress  evolution  equation  in  B  (cf.  Moran  et  al.  (1990);  Simo 
(1998b)).  Recall  the  stress  mappings  in  equations  86-88,  which  when  we  take  the  material 
time  derivative,  leads  to  the  following  equations 
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To  complete  these  equations,  we  need  the  material  time  derivative  of  the  stresses  in  the 
intermediate  configuration,  S^i,  an<l  Mrlm,  as  well  as  the  mapping  of  the  plastic 
evolutions  equations  to  the  current  configuration.  Given  the  constitutive  equation  for  S^l 
in  equation  116,  we  can  write 
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where  E\^  =  Ce^^/2  and  We  can  show  that 
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where  de  is  the  symmetric  clastic  deformation  rate  in  B,  and  £e  a  mixed  micro-macro 
clastic  velocity  gradient  in  B.  Then  we  can  write 
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where  the  left  elastic  Cauchy-Green  tensor  is  b\j  =  and  a  mixed  elastic  deformation 

tensor  ?/T  =  F^x^fq-  It  is  then  possible  to  write  the  material  time  derivative  of  the  stress 
difference  map  as 


KkVrl  -  Skl)f;l  =  ^4^  +  vbiftXfl 


For  the  couple  stress,  the  material  time  derivative  is  written  as 
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We  can  rewrite  the  gradient  of  elastic  micro  deformation  as 


such  that 


r %PQ  =  KNinPqxePpFeqQ  ,  Ienpq  =  x%xenA,q 


re-  -  - 

1  NPQ 


Fe-  y  ve-Fe- 

nN  1  nPq  X pPA  qQ 
°e  def 

7  =  Xe  +  Pe  rf  i  ,y 

'  nva  Inna  1  an  lava  1  n 


nPq  /  npq  1  an  /  apq  1  /  npa  aq  '  i  naq  aP 


+  7e  Ve 

1  /  naa  a 


Then,  the  couple  stress  material  time  derivative  map  becomes 
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where  bx’emp  =  X^AiXpM-  In  summary,  we  have  the  stress  evolution  equations  in  B  as 
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J~e  i(A  +  T)(beijdij)bekl  +  2(h  +  o)be kidljtf ji  + 

+Kbeki£eij'ipejl  +  (160) 

je 

Skl  &kl  TZ^kl  &kl)  T  1'  ki($il  0"iz)  T  &ki)^  li  T 

■J 

j;( Kbeu£^ij}ejk  +  ube  ki£ei:j^e  jt)  (161) 

je 

ril him  -j—Wlklm  “1“  ^  ki^Ulm  ^Y^kim^  li  ^Tlkli^  mi 

4  h  (vk,rnmrqp + r,mrvb\,)  +  r2  {b%,^mrwm + rkmb\vnr) 

+T3bekirqmrnp  +  n ^elmbekn^eqp  +  r5  {^ekmbeinVqp  +  VlmVkptfnq) 
+r6^ekmVipbenq  +  r7bekn^elp^eqm  +  r8  (V kpbe iqV nm  +  bekqbelnbx’emp) 

+T9beknbeiqbx’emp  +  TW^ekpbein^eqrn  +  Tii&%-0e^enm]  7np(?  (162) 

where  objective  elastic  stress  rates  are  defined  as 

°ki  =f  °ki  -  £ekicrii  -  Vkdu  +  deuaki  (163) 

7 - ^ - r  d=  {ski-dki) -tki(sii-aii)-(ski-aki)tu  +  deii(ski-aki)  (164) 

TYl him  dlkim  ^k i'nT'ilm  t^Inmd-n  ^'kli^rni  "b  did ^klrn.  (165) 

where  Je/Je  =  df{.  The  stress  rates  on  a  and  (s  —  a)  are  recognized  as  the  elastic  Truesdell 

stress  rates  (Holzapfel,  2000),  whereas  the  stress  rate  on  the  higher  order  stress  is  new,  and 

□ 

can  similarly  be  defined  as  an  elastic  Truesdell  higher  order  stress  rate.  To  show  m  is 
objective,  consider  a  rigid  body  motion  (Holzapfel,  2000),  with  translation  c  and  rotation 
Q  (orthogonal:  QQ2  =  1)  on  the  current  configuration  B,  resulting  in  B+,  such  that 

x/+  =  c  +  Q(x  +  £)  =  c  +  Qx'  (166) 

Recall  the  definition  of  the  higher  order  stress  through  the  area-averaging,  but  now  on  the 
translated  and  rotated  configuration  B+  as 
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def 


mklrnnkda  =  /  °kl  Zmnkda 

J  da 


where  o'^  =  Qka<j'abQib  ,  C  =  Qm£c  ,  n'£  =  Qkdn'd 

[  QuKiQM  me  tcQkdn'dda' 


J  da 
Qka 


'da 


O ab  £cndda  ^  QibQmcQ 


kd 


def  , 

—  'mabc'ft'dda 


QkaWlabcQlbQmcH'bd® 


(167) 


We  employ  the  standard  results  (Holzapfel,  2000) 


ifl 

T  Qki^ijQlj 

(168) 

Kt 

T  Qki^ijQlj 

(169) 

dtt 

=  dii 

(170) 

where  f =  QkaQia ■  We  substitute  into  the  expression  for  m  ,  with  after  some  tensor 
algebra,  we  can  show  that  m  is  objective: 


mklm  =  ™tlm  ~  ZeHmtlm  ~  “  mtliVmi  +  dtt mtlm 

Qka(rilklm  ^ki^ilm  ^Tlkim^li  ^'kli^rni  d,n'^'klm)QlbQmc 

Qka  TTlkl  m  QibQ  me 
q.e.d. 


(171) 


Now,  to  map  the  plastic  evolution  equations,  we  start  with  the  macroscale  plasticity.  The 
yield  and  plastic  potential  functions  in  B  become 


F(a,  c ) 


G(a,  c ) 


Je||deva||e  -  Je  (A*c  -  B^p)  <  0 
||devcr||e  =  (deva^&^V-^deva,™) 


devcTy  =  (T-,j 


r,  @kk  )  djj  ,  P  (jk k 


Je\\deva\\e  -  Je  (A^c-B^p) 


(172) 


(173) 
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where  be ^  =  Fe  d.Fe  d..  The  map  of  the  plastic  velocity  gradient  and  strain-like  ISV 
become 


4 


z 

c 


7 


3G 
dcr  id 
dG 


=  r 


if-1 
0  ka 


devcra6  x 


daH 

.  dG  ^ . 

—7—  = 

oc 

HZ 


1 1  dev  a  | 


+  -BHu 


(174) 


(175) 

(176) 


Next,  for  microscale  plasticity,  the  yield  and  plastic  potential  functions  in  B  become 


Fx(s  —  a,  cx) 

=  Je||dev(s  -  a)||e  -  Je  {Ax^cx  -  EX’V)  <  0 

1 

Px  =  -{Ski  ~  aki)5kt 

(177) 

Gx(s  —  a,  cx) 

=  Je\\dev(s  -  a)\\e  -  Je  (klx’V  -  5X’V) 

(178) 

The  map  of  the  plastic  microgyration  tensor  and  strain-like  ISV  become 


4  =  7X- 


dG x 


d(ski  -  akt ) 
dGx 


(179) 


=  Je 


Zx 


c '  = 


d{ski  -  aki) 
p)Gx 

-yx  =  Ax’^Je 
dcx 
HXZX 


( 6e-i  dev ( Sab  (Tab)  je-l  , 

\  ka  ||dev(s  —  cr)|p  bl  '  3  *kl 


(180) 

(181) 


For  microscale  gradient  plasticity,  the  yield  and  plastic  potential  functions  in  B  become 


FVx( m,  cVx) 


GVx(m,cVx) 


Je||devm||x  -  Je  (AVx^||cVx||x  -  SVx^||pVx||x)  <  0  (182) 

||cVx||x  =  \J &x^nCnX 

||devm||x  d=  (d evmijk)be^be~^bx’e^ (devmmnp) 

IIpVxIIx  =  ^Pmxbx^npnx  ,  p?x  =  ^ mkkm 

Je||devm||x  -  Je  (^Vx^||cVx||x  -  SVx^||pVx||x)  (183) 
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The  map  of  the  gradient  plastic  microgyration  tensor  and  strain-like  ISV  become 


ylm,k 


Vx 

cl 


(7Vx) 


dG Vx 

dm 


klm 


<9GVx 

-(7Vx) 


/  dev?nabc  i  e_i  1  /jVx.hx  Pq *  ftx,e-i 

V||devm||x°  mc+3  ^  ||pVx||x°  ™ 


Vx 


X)/g'vX  r  - 

=  yfVx.V-^Vxye  |  6X,e-l 

||cVxHx 


(184) 


(185) 

(186) 


Remark  9.  It  is  reassuring  to  see  that  the  left  hand  sides  of  the  plastic  deformation 
evolution  equations  174,  179,  184  simplify  considerably  in  the  current  configuration  from 
the  forms  in  the  intermediate  configuration  133,  139,  144.  They  are  further  simplified  when 
assuming  small  elastic  deformations. 

Remark  10.  Solving  numerically  for  the  increment  of  ufm  k  leads  to  the  solution  of  the 
increment  of  and,  in  turn,  the  evolution  of  the  couple  stress  over  time.  We  see 
this  by  taking  the  spatial  derivative  of  the  microgyration  tensor  as 


^  lm,k  ^ lm,k 

(187) 

ulm,k  =  \XlKXKm] 

V  lm,k  —  llmk  ~  ^  lefl  amk  A  v  ami  lak 

(188) 

Remark  11.  With  the  definition  of  the  plastic  evolution  equations,  plastic  deformation 
can  be  calculated,  and  elastic  deformation  updated  to  calculate  the  stresses.  A  FE 
implementations  solve  these  equations.  We  take  advantage  of  the  small  deformation 
elasticity,  such  that  elastic  deformation  in  the  current  configuration,  be^\  beij,  etc.,  can  be 
approximated  by  the  second-order  unit  tensor  8ir 

To  illustrate  the  implementation,  we  apply  assumptions  for  small  elastic  deformation  in  the 
next  section. 

Small  elastic  deformation  and  Cartesian  coordinates  for  current  configuration 

B:  Assuming  small  elastic  deformation,  the  tensors  be\3x ,  lfl31  i/j ei3 ,  bx,eJx  ~  di3  and  Je  ~  1, 
when  multiplied  by  another  variable  that  is  not  8i3  or  1,  such  as  ^^Vlevcr^6^1  «  d ever*,;. 
Also,  the  rate  of  elastic  volumetric  deformation  is  Je /  Je  =  d^.  With  these  approximations, 
in  summary,  we  have  the  stress  evolution  equations  in  B  as 
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Ski  —  &kl  — 


'Wl'klm 


~(^ii)aki  +  ^kiaii  +  aki^u  +  (A  +  T)(d%)5ki  +  2(/i  +  a)dekl 

+l(£ii)^ki  +  K£li  +  u£ik  (189) 

~(deii)(ski  -  Vki)  +  tki(su  -  ou)  +  (ski  -  aki)tu  +  n£elk  +  U£ekl  (190) 

oe 

(d'iijTYlkim  H“  ^ki^ilrn  “1“  ^'kli^rni  ^ klmnpq  ^  npq 

C-klmnpq  Tl  {^kl^nm^qp  “1“  $lm$np$kq)  ^~2  {^kl^mp^nq  $km$lq$np) 
^kl^qm^np  “1“  ^~4 ^Im^kn^qp  “1“  ^~5  ^km^ln^qp  “I-  &lm$kp$nq) 

H- ^~6^km^lp^nq  H“  ^~7^kn^lp^qm  H“  7g  ($kp$lq$nm  H“  &kq$ln$mp) 

~^~^~9^kn^lq^mp  7*10 $kp&  ln$qm  ^”ll $kq$lp$nm  (191) 


The  yield  and  plastic  potential  functions  in  £?  become 


F(a,  c )  =  || deva ||  -  (A^c  -  B^p)  <  0 
||devcr||  =  \J (devrr,,  )(devrr,, ) 

dev<T„  =  Oij  -  (1  a**)  ki,  P=  \atk 
G(a,  c)  =  ||dev£r||  -  (A*c-B*p) 


(192) 


(193) 


The  map  of  the  plastic  velocity  gradient  and  strain-like  ISV  become 


dG  _  devcrkl  1 

o —  —  Tib - if  +  ?B^dki 

dcrki  ||devcr||  3 

■dG 

-7—  =  A*  7 


(194) 


(195) 

(196) 


Notice  that  £p  =  ^{d G/da)T .  Next,  for  microscale  plasticity,  the  yield  and  plastic  potential 
functions  in  B  become 


Fx( s  —  a,  cx)  =  ||dev(s  —  cr)  ||  —  (Ax’^cx  —  Bx,^px)  <  0 
P  ^('SfcA:  & kk ) 

Gx(s  —  a,  cx)  =  ||dev(s  —  cr)||  —  (Ax,^cx  —  Bx,^px) 


(197) 


(198) 
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The  map  of  the  plastic  microgyration  tensor  and  strain-like  ISV  become 


4  = 


AyX 


dG x 


d(sM  -  aki) 
dGx  =  dev(ski  -  akl)  +  1^^ 
d(ski-aki)  ||dev(s  —  cr)||  3 


(199) 


kl 


ZX  = 


CX  — 


-syX. 


dGx 


dc* 
HXZX 


—  Ax^^x 


(200) 

(201) 


For  microscale  gradient  plasticity,  the  yield  and  plastic  potential  functions  in  B  become 


FVx(m,  ccVx) 


GVx(m,  cVx) 


||devm||  —  (AVx^||cVx||  —  5Vx,l?i||pVx||)  <  0 
||rVxi|  _  ,  /  vx  vx 

|| ^  ||  —  Gm  Gm 

||devm||  =  J  (devmijk)(devmijk) 
devrriijk  =  mijk  -  -Sifmaak 

||pVx||  =  \/p7nXPmX  ,  p7x  =  \rnkkm 

||devm||  -  (j4Vx>^||cv*||  -  SVx,^||pVx||) 


(202) 


(203) 


(204) 


The  map  of  the  gradient  plastic  microgyration  tensor  and  strain-like  ISV  become 


V; 


,p  _ 
lm,k 


=  uyx)  dG^x  dc^x  =  devm^  +  \b*x,i>8uJi 
dmkim'  dmkim  ||devm||  3 


Vx 

m 


p 


Vxl 


(205) 


..vx 


f)Gv* 

Zx  =  -fdVx')— —  =  AVxv/’('dVxV 
a  17  ’  dc7x  17  J||cVx|| 


(206) 


CJX  =  H^XZX  (207) 

Solving  numerically  for  the  increment  of  ufm  k  leads  to  the  solution  of  the  increment  of 

and,  in  turn,  the  evolution  of  the  couple  stress  mkim  over  time.  We  see  this  by  taking 
the  spatial  derivative  of  the  microgyration  tensor  as 
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^Zrr^fc 


(208) 


^ lm,k  T  ^ lm,k 

Vlm,k  =  [XIKXkIu]  tk 

^ lm,k  Xlmk  ^ laXamk  ~f"  ^amPtlak 


Boxes  1  and  2  provide  summaries  of  the  stress  and  plastic  evolution  equations,  respectively, 
in  symbolic  form  to  solve  numerically  in  time.  The  details  of  the  symbolic  equations  have 
been  provided  in  index  notation  in  this  section.  The  numerical  time  integration  scheme  is 
presented  in  section  2.3.4. 
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Box  1.  Summary  of  stress  evolution  equations  in  the  current  configuration  in  symbolic  notation. 


Box  2.  Summary  of  plastic  evolution  equations  in  the  current  configuration  in  symbolic  notation. 
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2.5  Numerical  Time  Integration 


The  constitutive  equations  in  section  2.3.3  are  integrated  numerically  in  time  following  a 
semi-implicit  scheme  (Moran  et  ah,  1990).  We  will  solve  for  plastic  multiplier  increments 
Ay  and  Ayx  in  a  coupled  fashion  (if  yielding  is  detected  at  both  scales;  see  Box  9),  and 
multiplier  AyVx  afterward  because  it  is  uncoupled.  It  is  uncoupled  because  of  the 
assumption  that  quadratic  terms  in  equations  110  and  111  were  ignored,  leading  to 
uncoupling  of  the  higher  order  stress  m  from  Cauchy  stress  o  and  microstress  s,  whereas  o 
and  s  remain  coupled  (thus  coupling  7  and  yx). 

We  assume  a  deformation-driven  time  integration  scheme  within  a  FE  program  solving  the 
isothermal  coupled  balance  of  linear  momentum  and  first  moment  of  momentum  equations 
573  and  574,  respectively,  such  that  deformation  gradient  Fn+i  and  microdeformation 
tensor  Xn+i  are  given  at  time  tn+ 1,  as  well  as  their  increments  AF„+1  =  Fn+1  —  F„  and 
Ayn+i  =  Xn+i  —  Xn-  We  assume  a  time  step  A t  =  tn+ 1  —  tn.  Boxes  3-8  provide  summaries 
of  the  semi-implicit  time  integration  of  the  stress  and  plastic  evolution  equations, 
respectively,  in  symbolic  form. 

oe 

To  obtain  ye  in  Box  1  through  7  ,  we  use  equation  208  such  that 


Vv  =  Vue  +  Vvp 

Vi/  =  V  [tor1] 

Vz/e  =  ye  —  z/eye  +  ueT  0  ye 

=*►  7e  =  z/e y6  -  iseT  ©  7e  +  Vi/  -  Vz/P  (218) 


Recall  equation  158  which  gives  the  equation  for  the  objective  rate  of  ye  as 


Y  +  £eT7e  +  ^  +  rf  Q  ^  (219) 

which  appears  in  equation  211  in  Box  1,  and  in  Box  4  for  the  numerical  integration.  For 
V(x„+ 1  —  Xn)  and  Vyn+i  in  Box  4,  because  y  is  a  nodal  DOF  in  a  FE  solution  and  thus 
interpolated  in  a  standard  fashion,  its  spatial  gradient  can  be  calculated. 

Box  9  summarizes  the  algorithm  for  solving  the  plastic  multipliers  from  evaluating  the 
yield  functions  at  time  tn+\.  It  involves  multiple  plastic  yield  checks,  such  that  macro 
and/or  micro  plasticity  could  be  enabled,  and/or  microgradient  plasticity.  Because  the 
macro  and  micro  plasticity  yield  functions  F  and  Fx,  respectively,  are  decoupled  from  the 
microgradient  plastic  multiplier  7Vx,  we  solve  first  for  the  micro  and  macroplastic 
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multipliers,  as  indicated  by  (I)  in  Box  9,  and  then  for  the  microgradient  plastic  multiplier 
in  (II)  afterward.  Once  the  plastic  multipliers  are  calculated,  the  stresses  and  ISVs  can  be 
updated  as  indicated  in  Boxes  5-8. 

This  micromorphic  plasticity  model  numerical  integration  scheme  fit  nicely  into  a  coupled 
Lagrangian  FEformulation  and  implementation  of  the  balance  of  linear  momentum  and 
first  moment  of  momentum.  Such  work  is  ongoing. 
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Box  3.  Summary  of  semi-implicit  time  integration  of  Cauchy  stress  a  and 

microstress-Cauchy-stress  difference  (s  —  a)  evolution  equations.  (»)tr  implies  the  trial  value,  in 
this  case  the  trial  stress.  Results  of  the  semi-inrplicit  time  integration  of  the  plastic  evolution 
equations  in  Box  5  are  included  here. 


On+i  =  (1  -  tr(Aid®+1))crn  +  (A Uen+1)an  +  crn(AUen+1)T  + 

(A  +  r)tr(Atd£+1)l  +  2  (n  +  cr)(Afd,em)  +  r/tr(Ate^+1)l 

+«(A'  t£en+1)  +  u{A  teen+1)T  (220) 

(s  -  cr)n+i  =  (1  -  tr(Afd®+1))(s  -  a)n  +  (A^®+1)(s  -  a)n  +  (s  -  cr)n(At^+1)T  + 

K(At£en+l)T  +  v{Ateen+l)  (221) 


AtC+i  =  AUn+1  -  A Upn+1  (222) 

Atf'n+i  =  (AF„4i)Fn+1 

A<+1  =  (A,„+1)(r*y  .  + 

tr(Atd®+1)  =  tr(At£n+i)  -  tr(Af^)+1)  =  tr(At£n+i)  -  Bi\ Aq^.+i) 


A*4+i  =  Ati/*+1  +  At£el+1 

AtVn+i  =  Att'n+i  -  Afz/+1 

Atvn+l  =  {AXn+l)Xn+l 


dev(s  —  cr)tr 


+  -Bx’H 

3 


AK+1  =  (A7nV)(^f  ,  ^=||dev(s_a)1 

tr(Ate*+1)  =  tr(Aten+i)  -  tr(Ate£+1)  =  tr(Afen+i)  -  Bx^( Aq*+1) 
tr(Aten+i)  =  tr(Afz/„+i)  +  tr(At^n+i) 


(223) 
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Box  4.  Summary  of  semi-implicit  time  integration  of  higher-order  stress  m  evolution  equation. 
Results  of  the  semi-implicit  time  integration  of  the  plastic  evolution  equations  in  Box  5  are 

included  here. 


mn+i  =  (1  -  tr(Afd,em))m„  +  (AUen+l)mn  +  mn  ©  (A ttn+1)T  +  mn(A tv*+1)T  + 

ci(A  tln+1)  (224) 


oe 

At  7„_ 


At7n+1  +  (AtC+l)T7n  +  7n(AtC+l)  +  In  ©  (At<+i) 

At7n+1  =  (A^n+i)7n  -  {^tuen+l)T  ©  7n  +  AtVj/n+i  -  AtVz/£+1 
Atv t'n+l  =  V(xn+ 1  -  Xn)Xn+l  ~  (Xn+ 1  ~  Xn)Xn+l(VXn+l)Xn+l 

AtV^+1  =  (A7**)(f™r 


xVx,tr  _ 


devmtr 

||devmtr| 


Box  5.  Summary  of  semi-implicit  time  integration  of  plastic  evolution  equations  in  the  current 

configuration. 
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Box  6.  Elastic-predictor-plastic-corrector  form  of  semi-implicit  time  integration  of  stress  a 
evolution  equation  in  the  current  configuration. 


un+1  =  atr  -  (A7n+1)D*>’tr  -  (A7£+1)D*’*>’tr  (238) 

cjtr  =  (1  -  tr(AUn+i))an  +  {Atln+1)an  +  an{AUn+l)T  +  (A  +  r)tr(At4+i)l  + 

2 {/i  +  cr)sym(Ai4+i)  +  77  [tr(Atz/n+1)  +  tr(AUn+i)]  1 

-\-n  \Atvn+i  +  (At£n+i)T]  +  v  [(Atvn+i)T  +  Af4+i]  (239) 

DP,tr  =  _Bfan  +  +  +  (A  +  T)^l  +  2(M  +  a)sym(rtr)  +  7/5^1 

+7cftr  +  i/(ftr)T  (240) 

Dx’p’tr  =  rjB^l  +  Mrx’tr)T  +  z/rx’tr  (241) 


Box  7.  Elastic-predictor-plastic-corrector  form  of  semi-implicit,  time  integration  of 
microstress-Cauchy-stress  difference  (s  —  a)  evolution  equation  in  the  current  configuration. 


(s  —  o  )n+ 1 
(s  -  u)tr 


EP.tr 

E^’P’tr 


(s  -  u)tr  -  (A7n+1)EP’tr  -  (A7x+1)Ex’P-tr  (242) 

(1  -  tr(A*4+i))(s  -  u)n  +  (At4+1)(s  -  a)n  +  (s  -  a)n(AUn+l)T 
+7c  [(Ati7n_|_i)r  +  Al4+i]  +  v  [Aii^i+i  +  ( At£njri)T ]  (243) 

-Bi\ s  -  a)n  +  (rtr)T(s  -  a)n  +  (s  -  a)n ftr  +  K(ftr)r  +  z/rtr  (244) 
Kf.x.tr  +  ^(rX.tr  f  (245) 


Box  8.  Elastic-predictor-plastic-corrector  form  of  semi-implicit  time  integration  of  higher-order 


couple  stress  m  evolution  equation  in  the  current  configuration. 


mn+i 

=  mtr  -  (A7n+1)KP,tr  -  (A7*+1)K*’P’tr  -  (A7J_^1)KVx,P,tr 

(246) 

mtr 

=  (1  -  tr(At4+i))mn  +  (Atf„+i)mn  +  m,,  0  (At4+i)r  +  m„(Afi/n+i)T 

+c:  |(Ati/n+i)7®  -  (Atun+1)T  ©  in  -  7,'(  ®  (Atun+i)  +  AtVvn+i 
+  (At4+l)T7n  +  7n(^4+l)] 

(247) 

KP,tr 

=  B* mn  +  (rtr)Tmre  +  mn  ©  rtr  +  ci  [rtr7£  +  4(ftr)T] 

(248) 

K\-P’tr 

=  mnrx,tr  +  ci  [(rx,tl’)T7®  -  rx’tr  ©7^  +  7 n  ©  (rX’trf] 

(249) 

KVX.P>tr 

=  ci(fVx’tr)T 

(250) 

60 


Box  9.  Check  for  plastic  yielding  and  solve  for  plastic  multipliers. 


(I)  solve  for  macro  and  micro  plastic  multipliers  Ay  and  Ay*: 

Step  1.  Compute  trial  stresses  crtr,  (s  —  cr)tr,  and  trial  yield  functions  Ftr,  Fx,u' 

Step  2.  Consider  3  cases: 

(i)  If  Ftr  >  0  and  Fx,tr  >  0,  solve  for  Ayn+i  and  Ay*+1  using  Newton-Raphson  for 
coupled  equations: 

F(an+1,cn+i)  =  F(  Ayn+i,  Ay*+1)  =  0  (251) 

Fx((s-a)n+1,c*+1)  =  Fx(Ay„+i,  Ay*+1)  =  0  (252) 

(ii)  If  FtT  >  0  and  Fx,tY  <  0,  solve  for  Ay^.+i  with  Ay^ ,  1  =  0  using  Newton-Raphson: 

F(an+i,cn+1)  =  F(  Ay„+i,  Ay*+1  =  0)  =  0  (253) 

(iii)  If  FtT  <  0  and  Fx,tr  >  0,  solve  for  Ay*+1  with  Ayn+i  =  0  using  Newton-Raphson: 

F*((s-a)n+1,c*+1)  =  Fx(Ayn+1  =  0,Ay*+1)  =  0  (254) 

(II)  solve  for  microgradient  plastic  multiplier  AyVx,  given  Ay  and  Ayx: 

Step  1.  Compute  trial  stress  mtr  and  trial  yield  function  F^x,tr 

Step  2.  If  FVx,t r  >  0,  solve  for  AyJ^  using  Newton-Raphson: 

=  ^Vx(AyJ+\)  =  0  (255) 
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3.  Upscaling  from  Grain-Scale  to  Micromorphic  Elastoplasticity 


In  the  overlapping  domain,  the  continuum- scale  micromorphic  solution  can  be  calculated  as 
a  partly  homogenized  representation  of  the  grain-scale  solution  (figure  6)^.  This  is  useful 
for  fitting  micromorphic  material  parameters,  and  also  in  estimating  DNS  material 
parameters  when  converting  from  micromorphic  continuum  FE  mesh  to  DNS  in  a  future 
adaptive  scheme.  Thus,  a  micromorphic  continuum-scale  field  □micromorphic  jg  defined  as  a 
weighted  average  (over  volume  and  area)  of  the  corresponding  field  dP’3,111  at  the 
grain-scale,  which  is  written  as  follows: 


| — | micromorphic, vol  def  ^/j — jgrai  n  \. 


def 


JQavs 

| — |micromorphic,area^|  def  / 1— |grain  j^grain  \  def  ^ 

\  /a  pavg 


u(r,  e ,  tf)Dgraindu 


□gramngramda 


fpavg 


(256) 

(257) 


where  (•)v  denotes  the  volume-averaging  operator,  rA,avg  =  /Qavg  u(r,  9,  i))dv  the  weighted 
average  current  volume,  uj(r,6,'&)  denotes  the  kernel  function  (if  using  spherical  coordinates 
in  3-D  averaging),  f2avg  is  the  grain-scale  volume  averaging  domain,  (•)  denotes  the 
area-averaging  operator,  and  Tavg  is  the  grain-scale  area  averaging  domain.  These  averaging 
operators  are  mapped  back  to  the  reference  configuration,  such  that  the  domains  DgVg  and 
TgVg  are  fixed.  A  length  i  (approximate  diameter  of  f2avg  and  Tavg)  is  a  material  property 
and  is  directly  related  to  the  length  scale  used  in  the  micromorphic  constitutive  model. 

A  macroelement  material  point  (figures  5  and  6)  can  be  characterized  as  fully  overlapping, 
non-overlapping,  or  partly  overlapping  according  to  the  level  of  overlapping  between  the 
averaging  domain  Davg  and  the  full  grain-scale  DNS  region  f2grain.  Within  the 
fully-overlapping  averaging  domain,  the  Cauchy  stress  tensor  u|)aiu  and  vector  of  ISVs  ggrain 
at  the  grain-scale  are  projected  to  the  micromorphic  continuum-scale  using  the  averaging 
operators  (□gram^  a,nd  (□grain)a  to  define  the  unsymmetric  Cauchy  stress  cr^,  the 


symmetric  microstress  Ski,  and  the  higher  order  stress  rrikim'- 

<yank  s  r<“)o  (258) 

=  {»tr)v  (259) 

mtlmnt  S  (260) 

^  (26!) 


^Wewould  like  to  acknowledge  discussions  with  his  colleague  at  University  of  Colorado  Boulder  (UCB), 
Prof.  F.  Vernerey,  regarding  the  natural  built-in  homogenization  in  micromorphic  continuum  theories  of 
Eringen  (1999);  Eringen  and  Suhubi  (1964). 
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where  it  is  assumed  the  variables  on  the  left-hand  sides  are  micromorphic.  Kinematic 
coupling  and  energy  partitioning  determine  the  percent  contribution  of  grain-scale  DNS 
and  micromorphic  continuum  FE  to  the  balance  equations  in  the  overlapping  domain. 


Figure  6.  Two-dimensional  illustration  of  micromorphic  continuum  homogenization  of 
grain-scale  response  at  a  FE  Gauss  integration  point  X  in  the  overlap  region. 
vRVE  implies  a  RYE  if  needed  to  approximate  stress  from  a  DE  simulation  at 
a  particular  point  of  integration  in  Davg,  for  example  in  Christoffersen  et  al. 
(1981);  Rothenburg  and  Selvadurai  (1981). 
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4.  Coupled  Formulation 


We  consider  here  the  bridging-scale  decomposition  (Kadowaki  and  Liu,  2004,  Klein  and 
Zimmerman,  2006,  Wagner  and  Liu,  2003)  to  provide  proper  BC  constraints  on  a  DNS 
region  to  remove  fictitious  boundary  forces  and  wave  reflections. 

4.1  Kinematics 

The  kinematics  of  the  coupled  regions  are  given,  following  the  illustration  shown  in  figure  2. 
It  is  assumed  that  the  micromorphic  continuum-FE  mesh  covers  the  domain  of  the  problem 
in  which  the  bound  particulate  mechanics  not  significantly  dominant,  whereas  in  regions  of 
significant  grain-matrix  debonding  or  intra-granular  cracking  leading  to  a  macro-crack,  a 
grain-scale  mechanics  representation  is  used  (grain-FE  or  grain-DE-FE).  Following  some  of 
the  same  notation  presented  in  Kadowaki  and  Liu  (2004),  Wagner  and  Liu  (2003), 
grain-FE  displacements  in  the  system  in  the  current  configuration  B  are  defined  as 


Q  =  ha,  q/3,  •  •  • ,  q7]T,  p,  ■  ■  ■ ,  7  e  A  (262) 

where  qQ  is  the  displacement  vector  of  grain-FE  node  a,  and  A  is  the  set  of  all  grain-FE 
nodes.  The  micromorphic  continuum-FE  nodal  displacements  da  and 
microdisplacement-gradients  <f>d  (see  below  for  xh  —  1  +  (Eringen  1968))  are  written  as 


D  =  [da,db, . . .  ,dc,(j)d,(j)e, . . .  ,(pf}T  (263) 

a.b. . . c  e  Sf  ,  d,  e, . . . ,  /  G  M. 

where  da  is  the  displacement  vector  at  node  a,  <f>d  is  the  microdisplacement-gradient  matrix 
at  node  d,  J\f  is  the  set  of  all  nodes,  and  A4  is  the  set  of  FE  nodes  with 
microdisplacement-gradient  DOFs,  where  typically  A 4  C  A/".  In  order  to  satisfy  the  BCs  for 
both  regions,  the  motion  of  the  grain-FE  nodes  in  the  overlap  region  (referred  to  as 
“ghost”  grain-FE  nodes,  figure  2)  is  prescribed  by  the  micromorphic  continuum 
displacement  and  microdisplacement-gradient  fields,  and  written  as  Q  G  A,  while  the 
unprescribed  (or  free)  grain-FE  nodal  displacements  are  Q  G  A,  where  A  U  A  =  A  and 
AnA  =  0.  The  displacements  and  microdisplacement-gradients  of  continuum-FE  nodes 
overlaying  the  grain-FE  are  driven  by  the  grain-FE  motion  (also  through  the  averaging 
shown  in  figure  6)  and  written  as  D  G  A f,Ai,  while  the  unprescribed  (or  free)  nodal 
displacements  and  microdisplacement-gradients  are  D  G  A f,M,  where 
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In  general,  the  displacement  vector  of  a  grain-FE  node  a  can  be  represented  by  the  FE 
interpolation  of  the  continuum  macro-displacement  held  and 

microdisplacement-gradient  held  <&h  (where  y/(  =  1  +  <&h  (Eringen,  1968))  evaluated  at  the 
grain-FE  node  in  the  reference  conhguration  Xa,  such  that 

uh(XQ,  t)  =  Y,  N2(Xa)da(t)  ,  <$>h(bXa,  t)  =  Y,  N*(Xa)<i>b(t)  a  e  A  (264) 

a&J\f  beM 

where  1V“  are  the  shape  functions  associated  with  the  continuum  displacement  held  uh, 
and  N®  the  shape  functions  associated  with  the  continuum  microdisplacement-gradient 
held  <&h.  Recall  that  iV“  and  Nf  have  compact  support  and  thus  are  evaluated  only  for 
grain-FE  nodes  that  lie  within  a  micromorphic  continuum  element  containing  nodes  a  and 
b  in  its  domain.  The  displacement  of  a  microelement  (figure  3)  can  be  written  as 


u'(X,  3,  t)  =  x'(X,  E,  t)  —  X'(X,  S) 

=  x(X,t)  +  £(X,E,t)-X- 
=  x(X,  t)  —  X  +  £(X,  S,  t) - 

u(X,i)  x(X,t)  3 

=  u(X,  t)  +  [x(X,  t)  —  1]  H 

" - V - ' 

*(x,t) 

=  u(X,t)  +  $(X,t)E 


(265) 


where  we  used  the  dehnition  y  =  1  +  $  (Eringen,  1968)  to  put  the  form  of 
microdeformation  tensor  y  similar  to  the  deformation  gradient  F  =  1  +  <9u/<9X.  The 
prescribed  displacement  of  ghost  grain-FE  node  a  can  then  be  written  as 

q a(t)  =  (u')/l(XQ,  Ea,  t)  =  uh(Xa,  t )  +  $ft  (Xa,  t) Ea  a  G  A  (266) 

where  Ea  is  the  relative  position  of  grain-FE  node  a  from  a  micromorphic  continuum 
point.  The  choice  of  this  continuum  point  could  be  either  a  continuum-FE  node  or  Gauss 
integration  point.  This  will  be  investigated  in  the  multiscale  implementation.  The  influence 
of  the  microdisplacement-gradient  tensor  <&h  in  the  overlap  region  on  the  ghost  grain-FE 
nodal  displacement  could,  through  specific  micromorphic  viscoelastic  constitutive  relations 
for  ye,  act  to  “damp  out”  high-frequency  waves  propagating  through  the  fine  mesh 
grain-FE  region  to  the  overlap/coupling  region.  The  partitioning  of  potential  and  kinetic 
energies  between  grain-FE  and  micromorphic-continuum-FE  systems  within  the  overlap 
region  will  be  dependent  on  the  grain-FE  equations  of  the  bound  particulate  system  and 
the  micromorphic  continuum-FE  equations  of  the  continuum  system. 
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For  all  ghost  grain-FE  nodes,  the  interpolations  can  be  written  as 


Q  —  •  D  +  N^g  •  D  (267) 

where  N and  N gg  are  shape  function  matrices  containing  individual  nodal  shape 
functions  iV“  and  N® ,  but  for  now  these  matrices  are  left  general  to  increase  our  flexibility 
in  choosing  interpolation/projection  functions  (such  as  those  used  in  meshfree  methods). 
Overall,  the  grain-FE  displacements  may  be  written  as 


'  Q  ' 

NQs1 

D 

+ 

'  Q' 

.  Q  . 

Nqd  _ 

D 

0 

where  Q'  is  introduced  (Klein  and  Zimmerman,  2006)  as  the  error  (or  “fine-scale”  (Wagner 
and  Liu,  2003))  in  the  interpolation  of  the  free  grain-FE  displacements  Q,  whose  function 
space  is  not  rich  enough  to  represent  the  true  free  grain-FE  nodal  motion.  The  shape 
function  matrices  N  are,  in  general,  not  square  because  the  number  of  free  grain-FE  nodes 
are  not  the  same  as  free  micromorphic-FE  nodes  and  prescribed  nodes,  and  number  of 
ghost  grain-FE  nodes  not  the  same  as  prescribed  and  free  micromorphic-FE  nodes.  A 
scalar  measure  of  error  in  grain-FE  nodal  displacements  is  defined  as  (Klein  and 
Zimmerman,  2006) 


e  =  Q'  Q'  (269) 

which  may  be  minimized  with  respect  to  prescribed  continuum  micromorphic-FE  nodal 
DOFs  D  to  solve  for  D  in  terms  of  free  grain-FE  nodal  DOFs  and  micromorphic 
continuum  FE  nodal  DOFs  as 


6  =  mbbn5b<«  -  N<^D)  •  mbb  =  nJbnoB  (270) 

This  is  known  as  the  “discretized  L2  projection”  (Klein  and  Zimmerman,  2006)  of  the  free 
grain-FE  nodal  motion  Q  and  free  micromorphic-FE  nodal  DOFs  D  onto  the  prescribed 
micromorphic-FE  nodals  DOFs  D.  LIpon  substituting  equation  270  into  equation  267,  we 
may  write  the  prescribed  grain-FE  nodal  DOFs  Q  in  terms  of  free  grain-FE  nodal  Q  and 
micromorphic-FE  nodal  D  DOFs.  In  summary,  these  relations  are  written  as 


Q  —  Bg^Q  +  BgDD  (271) 

D  =  BggQ  +  BgDD  (272) 
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where 


BQQ 

—  N  —  B~ 

lyQD°DQ 

(273) 

bqd 

~  N QD  +  NggBgD 

(274) 

B5q 

=  mAn^ 

DD  QD 

(275) 

Bdd 

=  -M~hN^N  QD 

DD  QD 

(276) 

As  shown  in  figure  2,  for  a  FE  implementation  of  this  DOF  coupling,  we  expect  that  free 
grain-FE  nodal  DOFs  Q  will  not  fall  within  the  support  of  free  micromorphic  continuum 
FE  nodal  DOFs  D  such  that  it  can  be  assumed  that  N qd  =  0  and 


(277) 

(278) 

where 


Q  —  B$qQ  BQDB 
D  =  B^gQ 


—  Nq5B5q 

(279) 

=  nqu 

(280) 

bBq 

=  M~hNT~ 

DD  QD 

(281) 

BBd 

=  0 

(282) 

The  assumption  ^  0  would  be  valid  for  a  meshfree  projection  of  the  grain-FE  nodal 
motions  to  the  micromorphic-FE  nodal  DOFs,  as  in  Klein  and  Zimmerman  (2006),  where 
we  could  imagine  that  the  domain  of  influence  of  the  meshfree  projection  could  encompass 
a  free  grain-FE  node;  the  degree  of  encompassment  would  be  controlled  by  the  chosen 
support  size  of  the  meshfree  kernel  function.  The  choice  of  meshfree  projection  in  Klein 
and  Zimmerman  (2006)  was  not  necessarily  to  allow  Q  be  projected  to  D  (and  vice  versa), 
but  to  remove  the  computationally  costly  calculation  of  the  inverse  M~h  in  equations  271 
and  272.  Since  we  will  also  be  using  the  TAHOE  code  for  the  coupled  multiscale 
grain-FE-micromorphic-FE  implementation,  where  the  meshfree  projection  has  been 
implemented  for  atomistic-continuum  coupling  (Klein  and  Zimmerman,  2006),  we  will  also 
consider  the  meshfree  projection  in  the  future. 

FE  Balance  Equations: 
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Following  standard  FE  methods  to  formulate  the  nonlinear  dynamic  matrix  FE  equations, 
using  the  DOFs  defined  in  the  previous  section,  the  balance  of  linear  momentum  for  the 
grain-scale  FE  is 


MqQ  +  F/ArT’Q(Q)  =  Fext’q  (283) 

where  is  the  mass  matrix  (lumped  or  consistent  (Hughes,  1987)),  F1NT,®(Q)  is  the 
nonlinear  internal  force  vector,  and  FEXT'®  is  the  external  force  vector  (which  could  be  a 
function  of  Q,  but  here  such  dependence  is  not  shown). 

For  the  balance  of  linear  momentum  and  balance  of  first  moment  of  momentum  in  equation 
57,  the  weak  form  can  be  derived  following  the  method  of  weighted  residuals  in  Hughes 
(1987)  (details  not  shown),  the  Galerkin  form  expressed,  and  then  the  FE  matrix  equations 
written  in  coupled  form  as 


MdD  +  F/7VT’d(D)  =  Fext’d  (284) 

where  M19  is  the  mass  and  microinertia  matrix,  F/JVT,Z)(D)  the  nonlinear  internal  force 
vector,  Fext,d  the  external  force  vector,  and  D  =  [d0]T  is  the  generalized  DOF  vector  for 
the  coupled  micromorphic  FE  formulation. 

These  FE  equations  can  be  written  in  energy  form  to  make  the  partitioning  of  energy  in 
the  next  section  more  straightforward.  For  the  FE  matrix  form  of  balance  of  linear 
momentum  at  the  grain-scale,  we  have 


d  fdT®\ 

di  V^0"J 


,  su“  pBXT,Q 

dQ  dQ 


where  T ®  is  the  kinetic  energy  and  U® 


the  potential  energy,  such  that 


(285) 


Tq  =  ^QM^Q 

UQ( Q)  =  /  Q  FINT'Q(S)dS 

Jo 

Carrying  out  the  derivatives  in  equation  285,  and  using  the  Second  Fundamental  Theorem 
of  Calculus  for  dU® /dQ,  leads  to  equation  283.  Likewise,  for  the  coupled  micromorphic 
balance  equations,  we  have  the  energy  form 


d_ 

dt 


(286) 


dTD\  dTD  |  dU°  -pEXT,D 


<9D 


<9D  <9D 


where  T'0  is  the  kinetic  energy  and  UD  is  the  potential  energy,  such  that 


Td  =  -DMdD 

2 

Ud(D )  =  /  FINT,D  (S)dS 

Jo 


4.2  Partitioning  of  Energy 


We  assume  the  total  kinetic  and  potential  energy  of  the  coupled 
grain-FE-micromorphic-FE  system  may  be  written  as  the  sum  of  the  energies 


T(  Q,D)  =  Tq(Q,Q(Q,D))  +  Td(D.D(Q))  (287) 

U( Q,D)  =  Uq(Q,  Q(Q,  D))  +  Ud(D,  D(Q))  (288) 


where  we  have  indicated  the  functional  dependence  of  the  prescribed  grain-FE  nodal  DOFs 
and  micromorphic-FE  nodal  DOFs  solely  upon  the  free  grain-FE  nodal  DOFs  and 
micromorphic-FE  nodal  DOFs  Q  and  D.  respectively.  Lagrange’s  equations  may  then  be 
stated  as 


d_  JdT\  _dT_  dU_ 
dt  \<9Q/  <9Q  +  <9Q 

d_  fdT\  _dT_  dU_ 

dt  V^dJ  ~ 


-pEXT,Q 

jpEXT,D 


(289) 


which  lead  to  a  coupled  system  of  governing  equations  (linear  and  first  moment  of  mo¬ 
mentum)  for  the  coupled  grain-FE-micromorphic-FE  mechanics.  Details  of  the  derivatives, 
partitioning  coefficients,  and  numerical  examples  will  follow  in  a  future  report  and  journal 
articles. 
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5.  Conclusion 


5.1  Results 

The  schematic  for  a  concurrent  multiscale  computational  modeling  approach  for  simulating 
dynamic  fracture  in  bound  particulate  materials  was  presented  that  accounts  for 
grain-scale  micro-cracking  influences  on  macroscale  fracture.  Details  of  a  finite  strain 
micromorphic  pressure-sensitive  Drucker-Prager  elastoplastic  constitutive  model  were 
presented,  as  well  as  its  semi-implicit  numerical  integration.  The  approach  for  coupling 
grain-scale  FE  equations  to  the  macroscale  micromorphic  FE  equations  was  presented. 

5.2  Conclusions 

A  three-level  (macro,  micro,  and  microgradient)  micromorphic  pressure-sensitive  plasticity 
model  provides  additional  flexibility  in  coupling  with  grain-scale  mechanics  in  an 
overlapping  region  (figure  1.2)  for  attempting  to  account  for  influences  of  grain-scale 
micro-cracking  on  macroscale  fracture  nucleation  and  propagation  under  dynamic  loading 
of  bound  particulate  materials.  The  thorough  formulation  of  the  finite  strain  micromorphic 
elastoplastic  constitutive  equations  in  the  context  of  nonlinear  micromorphic  continuum 
mechanics  has  been  established,  allowing  the  multiscale  framework  to  stand  on  a  firm 
footing,  which  heretofore  was  not  presented  in  the  literature. 

5.3  Future  Work 

Future  work  will  involve  completing  the  FE  implementation  of  the  finite  strain 
micromorphic  pressure- sensitive  Drucker-Prager  plasticity  model,  and  coupling  it  via  an 
overlapping  region  to  the  grain-scale  FE  mesh  where  a  projectile  may  impact  a  bound 
particulate  materials  target  (figure  2). 
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Appendix  A.  Derivation  of  F' 


The  formulation  of  equation  5  is  presented  in  this  appendix,  using  direct  notation.  To 
start,  we  recognize  that 


<9x'  _  <9x'  dX 
^X7  ““  <9X  9X! 


fix'  d'x  „  dE 

ax=F  +  ax=  +  x0x 

and 


(A-l) 

(A-2) 


<9X  i  dE 
M7  “  1  ~~  dX’ 

It  is  possible  to  show  that  dE/dX'  ps  dE/dX,  starting  with 


(A-3) 


dE  dE  dX 

dX1 7  “  dX  dX' 


(A-4) 


which  using  equation  A-3  leads  to 


<9E  /  dE\  1  dE 

^X7  =  V1  +  dXJ  dX 


(A-5) 


Assuming  the  gradient  of  microstructural  internal  length  is  small,  ||<9H/<9X||  -C  1  for  the 
region  of  interest  where  the  micromorphic  continuum  model  is  used,  then 


(A-6) 


where  then 


dE  f  dE\  dE  dE 
dX7  =  \~dXjdX~dX 

The  expression  for  F'  then  results  as  in  equation  5. 


(A- 7) 
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Appendix  B.  Another  Set  of  Elastic  Deformation  Measures 


Here,  another  set  of  elastic  deformation  measures,  from  equation  1.5.11  in  (Eringen,  1999), 
are  considered  as 


/-06e 

UKL 


vc  -  v  -  T  c--  =  v-h- 
X-kKX-kL  ’  J-iCL  A  Lar aK 


Thus,  the  Helmholtz  free  energy  function  is  written  as 


(B-l) 


MCfl,  T |.J,  Z,t,  Z|,  ZJ  J,  9) 


and  the  constitutive  equations  for  stress  result  from  equations  82  -  84  as 


(B-2) 


’KL 


Mrlm  ~ 


D'V'e  ^  Bk  Lk 

U  1  KB 

(B-3) 

te_l  d{p^)-e_l 

AR  dCx/B  BL 

(B-4) 

0{p^)  e—1  ye— 1  _ 

A  IkJr  Lk 

(B-5) 

where  Y^”-1  =  Fe_1  KkXkA'  These  stress  equations  take  a  somewhat  simpler  form  than  in 
equations  91  -  93.  Thus,  it  becomes  a  choice  of  the  modeler  how  the  specific  constitutive 
form  of  the  elastic  part  of  the  Helmholtz  free  energy  function  is  written,  i.e.  in  terms  of 
equation  90  or  equation  B-2.  Eringen  (1999)  advocated  the  use  of  equation  B-2  for 
micromorphic  elasticity,  whereas  Suhubi  and  Eringen  (1964)  used  equation  90.  We  use 
equation  90. 
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Appendix  C.  Deformation  Measures 


It  was  mentioned  that  the  change  in  the  square  of  microelement  arc-lengths  ( ds ')2  —  ( dS ')2 
should  include  only  three  unique  elastic  deformation  measures  (the  two  sets  proposed  by 
Eringen  (1999)  and  considered  in  this  report  for  finite  strain  elastoplasticity).  Here,  we 
write  directly 


{ds')2  =  dx'dx'  =  dx'kdx'k 


(C-l) 


where 


dx'k  ~  ^kkdXK  +  xekK,L^KdXi  +  xekRxPKK^KdXL  +  (C-2) 

Then 


_  I  pe_  _  _  rie-l-ne _ 77  _ 

•B  '  L  DAK^  DM  mbl^a^b 


‘A^E 


{ds')2  =  [C%L  +  2sym(F^£) 

+2sym(y%ECefdTeCARxlEL) 

+  ^AD^AB^SE  XPdd,k  ^EE,L  ED^E 

+2sym(T^x|BjZ)5p  dXRdXL 

+2  [^RL  +  ^SL^BO^CAk^A 

+nLcrinBxpBD,Rzn\  dxRdz-L 

+  [^Ak^ab^Sl]  dZRdZL 


(C-3) 


and 


(, dS 'f  =  dXRdXR  +  2dXRdZR  +  dZRdZR  (C-4) 

It  can  be  seen  that  the  first  set  in  equation  89  appears  exclusively  as  elastic  deformation  in 
equation  C-3;  there  are  also  some  plastic  terms,  which  do  not  appear  in  equation  1.5.8  in 
Eringen  (1999).  Equation  C-3  could  likewise  be  expressed  in  terms  of  the  elastic  set  in 
equation  B-l.  But  one  or  the  other  set  is  unique,  as  outlined  by  Eringen  (1999)  for 
micromorphic  elasticity,  here  put  into  context  for  finite  strain  micromorphic  elastoplasticity. 
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